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FOREWORD 

This document presents the results of work performed by 
the Fluid Mechanics Section of the Aeromechanics Department 
of Lockheed's Huntsville Research & Engineering Center. This 
report is Volume I of a four -part final report, as required to 
fulfill Contract NAS7-761. This work was sponsored by the 
Liquid Propulsion Section of Jet Propulsion Laboratories, Mr. 
Wolfgang Simon, Technical Manager. 

This document constitutes Volume IV of a four-part final 
report. The other three volumes, printed separately, are: 

Volume 1— "Summary Volume — Method of Char- 
acteristics Plot Program," LMSC-HREC D162220-I 

Volume II — "User's Manual — Method of Character- 
istics Plot Program," LMSC-HREC D162220-II 

Volume III — "Solution of Non-Isoenergetic Super- 
sonic Flows by the Method of Characteristics," 
LMSC-HREC D 162220 -III. 

REVISION NOTICE 

Lockheed Missiles & Space Company, Inc., Huntsville 
Research & Engineering Center Technical Report HREC-7761-I 
(D162220-I), entitled "User's Manual — Variable O/F Ratio 
Method of Characteristics Program for Nozzle and Plume 
Analysis," dated June 1971 is revised as indicated below. 

Revision A changes the document number from HREC- 

7761-1, D162220-I to HREC-7761-4, D162220-IV. 

Revision A affects title page and pages ii, 2-3, 2-4, 2-6, 2-7 

3-2, 3-3, 3-4, 3-5, 3-7, 3-9, 3-10, 3-11,3-13, 3-15, 3-16, 3-17, 

3-22, 3-24, 3-25, 3-26, 3-56, 3-71, 4- 1, A- 15, A- 16, A- 17 
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Symbol 


SYMBOLS AND NOTATIONS 
Description 


M 

R 

R, r 
X, x 

y 

p>p 

T, t 

S, s 
V , v 
M 

, * 
A/A 

w 

A 

y 

6 

€ 

9 

M 

0 

p 

X 

* 

<P 

Subscripts 

00 

o 

B 


Mach number 
"gas constant" 
radial coordinate 
axial coordinate 
normal coordinate 
pressure 
temperature 
entropy 
velocity 
mass flow 
area ratio 
weight flow 
area 

isentropic exponent 

turning angle through shock 

shock angle 

flow angle 

Mach angle 

characteristic angle 

density 

equation modifier 
operation function 
point description information 

freestream conditions 
stagnation conditions 
boundary conditions 
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Section 1 
INTRODUCTION 

A precise knowledge of local flow properties in nozzles and exhaust 
plumes is necessary for performance, radiation, attenuation, heat transfer 
and impingement analyses. All of these analyses are dependent on an accu- 
rate knowledge of the environment. 

Lockheed Missiles & Space Company, Huntsville Research &. Engineering 
Center has developed, under the sponsorship of several governmental agencies, 
a two-dimensional or axisymmetric method-of-characteristics program. The 
program is applicable for problems involving supersonic flow of an inviscid, 
adiabatic reacting gas in thermal equilibrium. 

Areas of particular interest are modifications to the existing program 
to expand the capabilities and flexibility of the program. In particular the 
capability to handle O/F gradients within the flowfield has been incorporated. 

As a preliminary step in the plume analysis, an accurate knowledge of 
the nozzle flow field is required. Both nozzle and plume calculations are per- 
formed with the same program. This program is a versatile, user-oriented, 
analytical tool which is capable of producing all of the gas dynamic nozzle or 
plume data required for an impingement analysis. The user may choose a 
solution for a 

1. nozzle only, 

2. plume only, 

3. combination nozzle and plume, 

depending on the options and starting data selected. This program has been 
in use for several years and experience has led to a continual refinement of 
the calculational procedure. 
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This document was prepared to £acilitate operation and under standing 
of the program. Descriptions of the individual routines are presented. 
Answers to any questions pertaining to the operation of the program should 
be found in this document or in the program listing. Questions involving 
initial assumptions made in applying the general theory can be answered by 
referring to Volume III of this document or to the Appendix. 
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Section 2 

PROGRAM CAPABILITIES AND APPLICATIONS 


2. 1 CAPABILITIES 

The nozzle and plume computer program described in this document 
can be used to solve a wide variety of problems in real gas, supersonic, 
compressible flow. Capabilities were previously discussed in Ref. 1; how- 
ever, as improvements continue to be made to the basic program new 
capabilities evolve. Some of the more important, basic capabilities of the 
existing program are outlined below: 

• The gas may be ideal or real. If real, frozen or equilibrium 
assumptions can be made. Oxidizer/fuel gradients may be 
considered. 

• Two-dimensional or axisymmetric problem geometry can be 
used. 

• Both upper and lower boundaries can be solid or free. (A 
solid boundary can be approximated by either a conic or poly- 
nomial equation. ) 

• A nozzle wall may be curve fit with discrete points. 

• One compression corner on the upper wall can be calculated. 

(Any number may be considered if the problem is re-started 
each time. ) 

• The number of Prandtl -Meyer rays to be computed around 
expansion corner discontinuities may be input. 

• Any number of expansion corners can be considered on either 
the upper or lower wall. 

• Various methods for obtaining an initial start line are utilized. 

1. The program will calculate a one -dimensional start line 
anywhere in the nozzle. 

2. The program will calculate a start line at points within 
the nozzle necessary to conserve mass. 

3. Characteristic data can be input at points across the flow 
field within the nozzle or in the plume. 
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4. Any right running characteristic line can be used for a 
start line. (See page 2-10.) 

5. Any left running characteristic line can be used (may be 
in combination with a normal start line). (See page 2-10.) 

6. Any left running line may be input with a right running 
shock crossing it. (See page 3-15.) 

• Hypersonic or quiescent approach flow options may be used. 

• Exit to ambient pressure ratios from over expanded to highly 
under expanded are possible. 

• Viscous boundary layer approximations at the nozzle lip are 
available. (See page 3-25 and Appendix.) 

• Displacement of the axis of symmetry from the center of flow 
(i. e. , the plug nozzle flow field) is possible. 

Reacting gas solutions have been facilitated by modifying the Chemical 
Equilibrium Composition Program (CEC) (Ref. 2) to provide binary tape and 
punched output of its equilibrium or frozen real gas calculations at any desired 
O/F ratio(s). The Method -of-Characteristics program has the capability for se- 
lecting the proper case from a large set of real gas properties cases stored on a 
master tape. Lockheed/Huntsville's master tape presently consists of approxi- 
mately 70 cases and is continually being expanded. The method of generating this 
master tape is outlined in the diagram on the following page. Cases stored 
are uniquely identified by some characteristic of the particular gas under 
consideration. For example, a LOX/LH^ system may be identified by the 
following: 


Gas Type Mixture Ratio Chamber Pressure 

02/H2 O/F = 1.5 - 8.0 PC = 546.0 

New cases of general interest may be added to the master tape; however, 
ad hoc cases should be prepared on a separate tape. Tape preparation 
sequence and c cmmunication with the Method-of-Characteristics program 
is diagrammed on the following page. 

Once the method-of-characteristics solution has been obtained, the 
output tape may be used to map the flowfield using the method of characteristics 
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Plot Program described in Volume III. The output tape may be used by the 
Method-of-Characteristics Radial Plot Program (Ref. 3) which determines 
the radial variations of flowfield properties across the nozzle and plume flow- 
fields at constant axial stations. The Plume Impingement Program (PLIMP) 
(Ref.4)may also be run todetermine the effects of the rocket exhaust plume on 
objects immersed in the plume. Sequencing and communication of auxiliary pro- 
grams with the Method-of-Characteristics program is on the preceding page. 

Two-dimensional or axisymmetric solutions are selected by simply load- 
ing a control word in the program input iata. This integer (0 or 1) is then 
multiplied by the term containing (l/r) in the governing differential equation. 

By appropriate description of the flow boundaries, it is possible to change 
from a solid to free boundary on either the upper or lower walls. Conversely, 
it is not possible to change from a free to a solid boundary on either wall. 

Compression corners are allowed only on the upper wall. Because the 
program is restricted to a single shock within the flow field, the number of 
compression corners is limited to one. Also, right-running shocks only can 
be handled, thus no provision exists for compression corners on the lower 
wall. Remembering that a mirror image solution is possible, the problem 
can be inverted for this type of solution. Note that the program will fortuitously 
handle two shocks if the first shock terminates by intersecting a lower free 
boundary before the second shock begins. 

By choosing an input option, it is possible for the program to automati- 
cally reflect a shock which intersects the horizontal axis. The method is 
illustrated in Section 2.2 and explained in the Appendix. 

When a shock intersects the horizontal axis and the automatic reflection 
option has not been flagged, the solution may be restarted by inverting the prob- 
lem and "regularly" reflecting the shock. A shock is "regularly" reflected by 
inverting the problem and using start line options to initiate the reflected 
shocks (Sections 2.2, 3.1.1 and 3.1.2). When a restart is required for a shock 
which has terminated at a "Mach Disc," a boundary equation simulating a cylin- 
drical "pipe" must be used to approximate the subsonic region (Section 2.2). 
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Problem solid boundaries can be described by a set of discrete points. The 
set of points is spline fit with individual polynomials for each pair of points,- 
consistent with the polynomial boundary specification option of the program. 
Specific usage instructions for this option are included in the input Section 
(3.1.1). 

Shock wave calculations may be discontinued when the shock strength 
decreases to an insignificant value. The shock will continue to be calculated 
in its usual manner until the change in stagnation pressure becomes less 
than a given input percentage. After the criterion is satisfied, the shock 
will no longer be iterated. If a value for the percentage change is not pro- 
vided in the input, the shock will be computed in its usual manner until the 
percentage change drops below 0.1. (See Section 3.1.1 for description of input 
for this option.) 

The program computes nozzle thrust by integrating the pressure dis- 
tribution along the nozzle wall, including the thrust increment between the 
last characteristic line inside the nozzle and the Prandtl-Meyer expansion 
at the nozzle lip. The vacuum specific impulse and nozzle thrust coefficient 
at each nozzle wall point is also calculated. 

A mesh control option exists whereby the maximum and minimum char- 
acteristic mesh size can be controlled by input to the program. Mesh control 
increases program accuracy in large vacuum plume calculations. A descrip- 
tion of mesh control is contained in the Appendix. 

The streamline at each characteristic data point is calculated and stored 
on the output tape. 

2.2 APPLICATIONS 

The most common applications of this program are in the areas of 
standard iocket nozzles and axisymmetric plumes. Many other complex flow 
fields can be treated, however, if the user is familiar with the flexibility of 
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the program and its options. Several problems are discussed in the remaining 
paragraphs of this section to illustrate the versatility of the program. 

Consider the boundary conditions given below. This problem was run 
two-dimensionally (although it could just as well have been asixymmetric) 
for the purpose of program demonstration. 



The problem above, while being restricted to one compression corner, 
could have had more expansion corners illustrating a more complex case. 


Another problem depicting the program's versatility is the following: 
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The restriction that the flow remain supersonic (inherent in the method- 
of-characteristics solution) must be observed through the flow field. 

External flow can be simulated by specifying the necessary stagnation 
conditions and inserting a two-dimensional or axisymmetric object in the 
flow field. 
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The method of manually reflecting a shock "regularly" or through a 
"Mach reflection" is outlined below. 




Restart (regular reflection) 


V Solid Upper Boundary 
at Axis of Symmetry 




Restart (Mach reflection) 
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When the automatic reflection option is selected and a shock terminates 
at the lower boundary the program will automatically reflect the shock as out- 
lined in the diagrams below. 


Example of Automatic Restart 
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The subsequent sketches illustrate the start line options available to the 
user to begin the method of characteristics solution. Lines a, b and c are 
normal start lines obtained by specifying the same axial location of the upper 
and lower limits of the start line. Lines d, e and f are right-running start 
lines obtained by specifying different locations of the upper and lower limits 
of the start line. Lines g, h and i are left-running start lines. All three 
types of start lines may be input or calculated by specifying the Mach number 
along the start line, the area ratio (A/A ) at the start line, or by using con- 
servation of mass and the Mach number at the upper wall. Lines a, d and g 
may not be calculated by area ratio since they initiate from the nozzle throat. 




The options available to manually input a start line for a reflected sho^k 
case are explained and illustrated in Section 3.1.2 (ICON(9)). 
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Section 3 
DISCUSSION 

The program consists of 66 active subroutines and functions which per- 
form all the required calculations in the nozzle and plume. This building 
block approach simplifies troubleshooting, operation and modifications. 

Each of the 66 routines is explained separately. In addition, a general flow 
chart (Fig. 3.1) , and a breakdown of the routines according to function are 
presented. 

The input instructions and description of the Input FORTRAN symbols 
are covered in the input section. The output interpretation is covered in the 
output section. Also included is a list of all the common block variables 
along with a definition or use of each variable. 

3. 1 INPUT INSTRUCTIONS 

This subsection contains a detailed description of the program input 
instructions as follows: 


• Detailed, input guide 

• Detailed description of the input FORTRAN symbols 
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3. 1. 1 Detailed Guide for Input Data to MOC-O/F Program 

(Section 3.1.2 contains a detailed description of each of the input variables) 

CARD NO. 1 Problem Title or Identification 


Format: 
Cols. 1-72 


CARD NO. 2 


12A6 

HEADER(I) Comment card, header information such 
as problem title may go on this card. It 
is printed at the top of each page of output. 

Run Control Card 


Format: 


1615 (right adjusted) 


Col 5 
Col 9 


Col 10 


Cols 14, 15 
Cols 18*20 


ICON(l) 

ICON(2) 


ICON(2) 


ICON(3) 
ICON (4) 


1 Read cards for gas properties 

2 Read tape 10 (A6) for gas properties 

0 Regular start line (other than right or 
left running) 

1 Right-running characteristic start line 

2 Left- running characteristics start line 

0 Straight start line M given 

1 Source start line A/ A* given 

2 Starting line input 

3 Starting line calculated by conservation of 
mass using a linear Macn number distributioi 

Number of starting line points (50 max) 

Number of upper boundary equations including 
free boundary equation (100 max) 


Option for ICON (4) when upper boundary is to be curve fit 


Cols 17-20 ICON(4) 1000 + number of discrete points specifying 

upper boundary + 2, where the 2 represents 
the nozzle throat equation and the free 
boundary equation 

Cols 23-25 ICON(5) Number of lower boundary equations 

(100 max) 


Option for ICON(5) when lower boundary is to be curve fit 

Cols 22-25 ICON(5) 1000 + number of discrete points specifying 

lower boundary + 2, where the 2 represents 
the nozzle throat equation and the free 
boundary equation. 

Col 27 ICON (6) 0 The start line will not be output on 

punched cards 

1 The start line will be punched on cards 
in the form of Card 9. (This start line 
can be used to manually restart an 
"automatic" shock reflection if the pro 
blem terminates incorrectly after the 
shock is reflected by the program.) 
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(Card No. 2 Continued) 
Col 28-29 ICON(6) 


Col 30 ICON(6) 


N The boundary equation at which it is 
desired to begin the characteristic 
solution to automatically restart a 
reflected shock. This must be an 
upper solid wall equation (see Section 
3.1.2) 

0 Only the flowfield data generated after 
a shock is automatically reflected will 
be saved on the output tape. 

1 The data generated before a shock is 
automatically reflected will be saved on 
the first file of the output tape and the 
flowfield data generated after the shock 
is reflected will be written on the second 
file of the flowfield tape. (This option re- 
quires the use of computer system routines 
which space past records on data tapes.) 


NOTE: ICON(6) is the controlling flag for the "automatic’' shock reflection. 
Cards 13, 14 and 15 are also required for "automatic" reflection. 


Col 35 ICON(7) 

Cols 37-38 ICON(8) 


Col 39 ICON(8) 

Col 4U ICON(8) 


Cols 41-43 ICON(9) 


Cols 44, 45 ICON(9) 

Col 50 ICON(IO) 

Cols 51-55 ICON(ll) 

Col 60 ICON (12) 


0 Two-dimensional solution 

1 Axisymmetric solution 

This option controls the type of output ob- 
tained after the problem has reached a free 
boundary equation. The same scheme is 
used in cols 37-38 as in 39 and 40. If 
blanks are used in Cols 37-38 then Cols 39 
and 40 will control the printout for the en- 
tire run. (Ref. Section 3.1.2). 

0 Full output 

1 Limited output (Boundary, shock, input and 
P-M points) 

1 One line output 

(R, X, M, THETA, S, Shock Angle) 

2 Two lines, above plus 
(Mach Angle, P, Density, T, V) 

3 Three lines, aboye plus 
(MWT, GAMMA, TO v , PO , S ) 

No. of left -running points up to and in- 
cluding upstream shock point. Used when 
ICON(2)>20 and shock crosses starting 
line. (See Section 3.1.2, p.3-14.) 

Number of regular start line points if 
ICON(2) >20 
0 Not presently used 

Case number (prints at top of each page) 

0 Calculate shock wave 

1 No rotation option 
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(Card No. 2 Continued) 

Col 65 ICON(13) 0 

1 


Cols 69-70 

ICON (14) 

N 



O 

Col 74 

ISTOP 

0 



1 

Cols 75 

ICON(15) 

0 



1 

Col 80 

ICON (16) 

0 


1 


Do not use the viscous option to set 
up the start line 

Use the viscous option to set up a start 
line at a nozzle exit plane. (See Card 12, 
Section 3.1.2 and Appendix)* 

Number of Prandtl-Meyer rays to be 
used in expansions 

Program will determine number of rays 
to be used. 

Single case is being run 
Multiple cases are being run 
Mesh control is not desired 
Mesh control is desired 
(Ref. Section 3.1.2 and Appendix) 

No printout +!(e 

Printout intermediate data 


Entire start line must be a regular start line at the exit plane. 

j)C 

This option causes intermediate calculations to be printed out during shock 
and pitot pressure calculations. This option was used during program check 
out and is not generally required except for debugging purposes. 
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CARD NO. 3 Describes physical boundaries of the flow 

field. For use when IC0N(4) and/or 
ICON(5) < 100. 

Format: II, 3X, II, 

5X, 6E10.6 


Col 1 

Col 1 
Col 1 


Col 1 


IWALLi 1, Conic equation u 

R=A* (SORT (B+C V X+D X 2)+E) 

(See Section 3.1.2 pg. 3-17 for an example 
and description) 

IWALL 2, Polynomial equation 

R=A*X*$ 4 + B*X**3+ C*X**2 + D*X + E 


IWALL 3, Free boundary equation 

P=PINF* (1+E*X)*(1+GAMMAINF*(MINF*SIN*(THETAB 
- THETAINF))**2 

(See Section 3.1.2 pg. 3- 17 for an example and description) 
IWALL 6, Free boundary equation of same form 

as IWALL =3, except oblique shock solution 
is used to generate plume boundary. Use 
number 6 for free stream approach flow in 
range of Mach 1.5 to 5.5. 


NOTE: The problem uses air with a molecular weight of 28.966 for these 
oblique shock calculations. 


Col 5 


Cols 11-20 
21-30 
31-40 
41-50 
51-60 
61-70 


ITRANS 0 No discontinuity follows this equation 

ITRANS 1 Expansion corner follows this equation 

ITRANS 2 Compression corner follows this equation 

WALLCO(I, 1, J) A(If IWALL = 1 or 2), PINF (If IWALL=3) (Psfa) 

WALLCO(I, 2, J) B(If IWALL = 1 or 2), GAMMAINF (If IWALL =3) 

WALLCOU, 3, J) C(If IWALL = 1 or 2), MINF (If IWALL=3) 

WALLCO(I, 4, J) D(If IWALL = 1 or 2), THETAINF (If IWALL=3) 

WALLCO(I, 5, J) E(If IWALL = 1 or 2), E(If IWALL=3) 

XMAX Maximum X value for which this equation 

applies. 


NOTE: The coefficients of each equation are contained on a single card. 

As many cards, i.e., equations, as necessary to describe the 
boundaries are input. The units of physical dimensions affect only 
the thrust calculations in which units of feet are assumed. Upper 
boundary information is given first. Program assumes that starting 
line is bounded by solid walls and that the equations are ordered 
with XMAX monotonically increasing. 
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CARD NO. 3 A 
Fo rmat: 

6E10.6 

Optioned, physical flowfield boundary con- 
dition for use when ICON(4) and/or ICON(5) 
> 1000. 

Cols 1-10 

RC 

Radius of curvature of nozzle throat 

11-20 

RT 

Radius of nozzle throat 

21-30 

THETA 

Nozzle throat divergence angle 

31-40 

XO 

Axial coordinate shift (Axial distance from 
the origin of the coordinate system to the 
throat) 

Cols 1-10 

XIN(I) 

Each pair of XIN(I) and YIN(I) represents 

11-20 

YIN(I) 

a discrete point for spline -fitting of the 

21-30 

XIN(I+1) 

nozzle solid boundary. As many of these 

31-40 

YIN(I+1) 

points as necessary can be input to a 

41-50 

XIN(I+2) 

maximum of 100. Points must be input 

51-60 

YIN(I+2) 

with XIN(I) monotonic ally increasing. 

NOTE: The first point describing the wall must be downstream of the 
throat equation* 

Cols 1-10 XIN(I+3) 

YIN (1+3) 
etc 

Format: 

• 

11, JX.Il, 
5X, 6E10.6 

— 

Col 1 

IWALJL 

3 Free boundary equation 

Col 5 

1TRANS 

0 No discontinuity follows this equation. 

Col 11-20 

WALLCO 

PINF 

21-30 

II 

GAMMAINF 

31-40 

II 

MINF 

41-50 

11 

THETAINF 

51-60 

11 

E 

61-70 

II 

XMAX 


NOTE: Optional boundary description, Card Type 3A, is set up specifically 
for a rocket nozzle. The first card describes the nozzle throat, 
which is followed by cards containing sets of discrete points, tl-.ree 
per card, describing the nozzle contour. The points will be auto- 
matically spline-fit to fcrm equations for the contour. It is assumed 
that a free boundary equation follows the last nozzle point. Card 
types 3 and 3A may be mixed, i.e., the upper boundary may be 
described by type 3A while the lower boundary is described by type 
3, or vice versa. Card type 3 and 3A may not be used simultaneously 
for a given boundary. 
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CARD NO. 4 


Format: 

4A6, 5X, A3, 
6X, 12, 3X, 12 

Cols 1-24 

ALPHA 

Cols 30-32 

UNITS 

Cols 39-40 
Cols 44-45 

IOF 

IS 

CARD NO. 5 


Format: 

E10.6 

Cols 1-10 

OFRAT 

CARD NO. 6 


Format: 

E10.6, 8X, 12 

Col 1-10 
Col 19,20 

STAB 
IV TAB 

CARD NO. 7 


Fo rmat: 

5E10.6 

Cols I -10 
Cols 11-20 

TAB(I, J, K, 1) 
TAB(I, J, K, 2) 

Cols 21-30 
Cols 31-40 
Cols 41-50 

TAB(I, J, K, 3) 
TAB(I, J, K, 4) 
TAB(I, J, K, 5) 


Gas Identification and Gas Property Input 
Control 


Gas name, identification for real gas 
properties on tape. May be any name 
when gas properties are input via cards. 

ENG English units are to be input (cards only) 
MKS Metric units (cards or tape) 

Number of o/F cuts, l.min., 10 max. 

Number of en*t*opy cuts (ignored for tape, 

1 for ideal ga», 2 max for real gas via cards). 

O/F value of each table (must be input even 
if constant O/F case is run). o/F tables must 
be input with O/F values monotonically increasin' 


Entropy value and number of velocity cuts. 
(Not input if ICON(l) = 2, i.e., gas proper- 
ties via tape; 0.0 if ideal.) 


Entropy value 

Number of Mach numbers for this entropy 
value 13 max, (1 if ideal gas) 

This card(s) gives the Mach number and 
associated gas properties at that Mach 
number and entropy. 


Mach number 

Gas constant (R) if UNITS = ENG, Molec- 
ular weight (MWT) If UNITS = MKS. 
GAMMA 

Temperature at this Mach number 
Pressure at this Mach number 


U nits 

Temperature 

Pressure 

Gas Constant (R) 


MKS 

°K 

Atmospheres 

(Molecular 

Weight) 


ENG 

°R 

psfa 

1545.2x 

32.2/Mwt 
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(Card 7 continued) 

NOTE: If two O/F tables were input with 2 entropy tables and 10 velocity 
"cuts" for each entropy the cards would be arranged as follows: 

Card 5 
Card 6 

Cards 7 (10 ea) 

Card 6 

Card 7 (10 ea) 

Card 5 
Card 6 

Card 7 (10 ea) 

Card 6 

Card 7 (10 ea) 

NOTE: Cards number 5, 6, and 7 are omitted if gas properties are input 
via tape. 
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CARD NO. 8 


This card specifies the necessary infor- 
mation for the starting line. (If ICON(2) 



/ 2, 12 or 22) (Ref. Section 3.1.2) 

Format: 

8E10.6 


Cols 1-10 

CORLiIP(l) 

Axial coordinate of upper limit of start 
line. 

Cols 11-20 

CORLiIP(2) 

Axial coordinate of lower limit of start 
line. 

Cols 21-30 

CORL,IP(3) 

Mach number or A/A* for start line. 

Cols 31-40 

CORLIP(4) 

Entropy level of start line. 

Cols 41-50 

STEP(5) 

Area of nozzle throat (units consistent 
with boundary equations) (Must be input 




for thrust coefficient calculation). This 
must also be input if the start line is 
generated by conservation of mass. 

Cols 51-60 

CORLIP(5) 

Constant O/F value 

Cols 61-70 

STEP(6) 

Minimum AP for discontinuing shock calc. 
If zero is used the program assumes 
.001 (Ref. Section 3.1.2) 

CARD NO. 9 


These cards are used to read in a known 
starting line (ICON(2) = 2, 12 or 22). 

Omit when Card no. 8 is used. As many 
of these cards are input as specified by 
ICON(3) on the run control card. Used 
only if ICON(2) = 2, 12, or 22. (Ref. Section 
3.1.2) 

Format: 

8E10.6 


Cols 1-10 

PSI(K, 1) 

Radial coordinate of this point 

Cols 11-20 

PSI(K, 2) 

Axial coordinate of this point 

Cols 21-30 

PSI(K,.3) 

Mach number of this Doint 

Cols 31-40 

PSI(K. 4) 

Flow angle of this point 

Cols 41-50 

PSI(K. 5) 

Entrcpy level of this point 

Cols 51-60 

PSI(K, 6) 

Shock angle of downstream shock point. 

Cols 51-60 

STEP(5) 

For last card only STEP(5) where STEP(5) 
is the area of nozzle throat (units consistent 
with boundary equations, must be input for 
thrust calculation). 

Cols 61-70 

STEP(6) 

For last card only STEP(6) where STEP(6) 
is the minimum A P for discontinuing shock 
calculation. 

Cols 71-80 

PSI(K, 8) 

Value of O/F at each starting line point. 


NOTE: To use the variable O/F capability, supply a value of O/F for each 
input point. Option 2 for ICON(2), must be selected for starting 
line input. For constant O/F operation supply one O/F table and 
set O/F constant on card 8. 
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CARD NO. 10 This card contains the necessary informa- 

tion to limit the calculations to those areas 
of interest. An unusual scheme is employed 
in order to make these limits efficient for 
the many problem orientations which are 
possible. Units should be consistent with 
those used on card 3. (Ref* Section 3.1.2) 


Format: 

8E10.6 


Cols 1-10 

CUTDAT(l) 

Radial coordinate defining upper cutoff. 

Cols 11-20 

CUTDAT(2) 

Axial coordinate defining upper cutoff. 

Cols 21-30 

CUT DAT (3) 

Angle cutoff line makes with horizontal. 

Cols 31-40 

CUT DAT (4) 

Radial coordinate defining downstream 
cutoff. 

Cols 41-50 

CUTDAT(5) 

Axial coordinate defining downstream 
cutoff. 

Cols 51-60 

CUTDAT(6) 

Angle cutoff line makes with horizontal. 

Cols 61-70 

STEP(3) 

Point Insert Criteria. If the change in 
axial distance between two characteristic 
lines along the axis exceeds Step (3) a new 
characteristic point will be inserted. If 
nothing is input the program will use 10,000. 

Cols 71-80 

TINF 

Free stream temperature (°R). Use only 
if IWALL=6 on card 3. 

CARD NO. 11 


Mesh control parameters, controls mesh 
size by inserting or deleting points. Use 
this only if ICON(15) on Card 2 is greater 
than zero. (Ref. Section 3.1.2 and Appendix) 

Format: 

6E10.6 


Cols 1-10 

DIV(l) 

Length control abs. AS insert 

Cols 11-20 

DIV(2) 

Length control abs. Af delete 

Cols 21-30 

DIV(3) 

Mach No. control abs. AM insert 

Cols 31-40 

DIV (4) 

Mach No. control abs. AM delete 

Cols 41-50 

DIV(5) 

Flow angle control abs. A0 insert 

Cols 51-60 

DIV (6) 

Flow angle control abs. A0 delete 

CARD NO. 12 


This card contains the input information 
for the viscous boundary layer option. Use 
this card only if ICON(13) on Card 2 is 
greater than zero. (Ref. Section 3.1.2 and 
Appendix) 

Format: 

11, 2X, 12, 
5X, 2E10.6 


Col 1 

N POWER 

Reciprocal of the exponent of the velocity 
profile in the boundary layer 

Cols 4, 5 

NBLPTS 

Number of boundary layer points specified 

Cols 11-20 

XL 

A characteristic length (usually nozzle wall 
length) 

Cols 21-30 

CU 

Conversion factor for mixed units of length. 
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Automatic Shock Reflection (Cards 13, 14, and 15) 

If the automatic shock reflection option is used a namelist must be used 
to change some of the original input data. It is possible to change any of the 
ICON variables of Card 2, cutoff data of Card 10, mesh control criteria of 
Card 11 and point insert criteria of Card 10. The cutoff data must be changed 
to indicate a reflected or inverted case. Angle cutoff parameters of Card 
10 are input in radians on the namelist. 


CARD NO. 13 (For reflected shock cases only) 


Col 2-7 $INPT2 Control card for namelist 


CARD NO. 14 


Format: Free Field 

Col 2-72 Changes to data consisting of for example: 

*Cutoff data array - CUTDAT(6) 

^Mesh control array - DIV(6) 

^Control array - IC0N(16) 

Point insert - STEP(3) 

A typical example of the namelist card is: 

CUTDAT(l) = 100., CUT DAT (2) = 0.0, CUTDAT(3) = 0.0 
CUT DAT (4) = -100., CUTDAT(5) = 100., CUTDAT(6) = 1.57 
DIV(l) =40., ICON (6) = 1000 

CARD NO. 15 


Col 2-5 $END End of data card. 


Cutoff data array is mandatory change to reset problem limits. All other 
namelist changes are optional. 
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3. 1. 2 Input FORTRAN Symbols 


Symbol 


Description 


Units (j appl) 


HEADER (I) 


ICON(l) 


ICON(2) 


ICON(3) 


This array contains the problem descrip- N/A 

tion which is written at the top of every 
page of printed output. This header in- 
put with Card 1 of the input data. 

Used to tell the program if the gas proper- N/A 
ties are read in from cards (ICON(l) =1) 
or read from tape 10 (A6)(ICON(l) = 2) 
which was generated by the NASA Lewis 
Thermochem Program and Tapgen program. 

After the gas properties have been written 
on the output tape ICON(l) is used by sub- 
routines TABLE and FABLE for communi- 
cation. 

ICON(2) is the type of starting line the N/A 

program is to calculate or read in from 
cards. ICON(2) consists of 2 digits. 

The first digit is a 0, 1 or 2 for a normal 
start line, a right-running start line or 
a left -running start line respectively. 

The second digit is a 0, 1, 2 or 3 for a 
straight start line Mach No. given, source 
start line A/A* given, starting line input 
or starting line calculated by conservation 
of mass respectively. The second digit is 
used in subroutine PLUMIN for a computed 
go to statement which calculates or reads 
in the starting line depending on the type 
desired. The first digit is used in 
PHASE 1 to initialize the starting line 
into a usable array for beginning the 
method-of-characteristics solution. 

The total number of starting line points N/A 

input or desired to be calculated by the 
program. Used as a reference for 
determining when the starting line 
points are used up. 
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Symbol 

ICON(4) 


ICON(5) 

ICON(6) 


Description 

ICON (4) is a dual purpose input para- 
meter. When it is not desired to 
input an upper boundary with points , 
ICON(4) is the number of upper boundary 
equations and ICON(4) is used to read in 
the desired number of upper boundary 
equations. When it is desired to curve- 
fit the upper boundary ICON(4) is 1000 
plus the number of discrete points speci- 
fying the upper boundary plus 2 where 2 
represents the nozzle throat equation and 
the free boundary equation. For the curve 
fit case ICON(4) is used to read in the 
information on the nozzle throat equation t 
the points representing the wall and the 
free boundary equation. 

Same as ICON(4) except ICON(5) applies 
to the lower boundary equations. 

This parameter supplies the program 
information necessary to automatically 
reflect a shock which intercepts the 
horizontal axis. ICON(6 l is a four-digit 
number. The first digit is a 0 if it is not 
desired to have the start line punched and 
a 1 if it is desired that the start line 
be punched. The second and third digits 
are the number of the upper boundary 
equation which the program will use in 
restarting the problem. The program 
uses the last point on this upper boundary 
as the first point on the starting line. The 
fourth digit is either a one or a zero. 

If a one is used, an end of file mark is 
placed after the MOC data already on 
the output tape. If a zero is used, the 
MOC output created after the shock is 
reflected is written over the original 
output which was created before the re- 
flection. This option requires system 
routines to space past end-of-file marks 
on flowfield tape which are not a part of 
the MOC program and must be a part of 
the individual computer system. 
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Symbol 

ICON(7) 


ICON(8) 


ICON(9) 


Description Units 

The type of flowfield solution which is N/A 

to be performed. A (0) for a two- 
dimensional solution and a (1) for an 
axisymmetric solution. ICON(7) is 
used in the mass flow calculations and 
in the axisymmetric term of the com- 
patibility equation used for the solution 
of each characteristic poinc. 

lCON(8) indicates the type of printed out- N/A 

put desired by the user. ICON(8) con- 
sists of two pairs of digits. The first 
pair controls the type of printout after 
the program begins using a free boundary 
equation for problem limits. The second 
pair indicates the type of output desired 
before the program begins using a free 
boundary equation or throughout the 
entire plume if nothing is input for the 
first two digits. The first digit of each 
pair controls the output at interior points 
on a characteristic line. If a zero is 
used the program will print out informa- 
tion at every characteristic point. If a 
one is used, printout will occur only at 
the boundary and shock points. The 
second digit is a 1, 2 or 3 for one, two 
or three line output at each point. (See 
output example. ) 

ICON(9) is a dual purpose control for N/A 

starting the characteristic solution of a 

reflected shock case. The first three 

digits are used for inputting the start 

line for a reflected shock case in which a 

left-running characteristic is crossed by 

the shock. These digits represent the 

total number of left-running points up to 

and including the upstream shock point. 

The last two digits are numbers of regular 
start line points. Two examples of this 
input are illustrated on the following page. 
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Upstream 
Shock Point' 


V 


Downstream 
Shock Point 



Regular Start 
Line Points 


ICON(9)=00107 

ICON(2)=22 

ICON(3)=8 


Upstream 
Shock Point 


ICON (9) =00 704 

ICON(2)=22 

ICON(3)=12 


Left Running 
Points 




Downstream 
Shock Point 


Regular Start 
Line Points 


Symbol 


Description 


Units 


ICON(IO) 


ICON(ll) 


ICON(12) 


A flag which tells the program if a N/A 

radiance tape is desired. A zero 
indicates no radiance tape and a 1, 
one tape, a 2, two tapes. The sub- 
routines for generating the radiance 
tapes are called through ICON (10) but 
the routines are presently dummied 
out. (See Capabilities section. ) 

The case number for the problem N/A 

under consideration. This number 
is printed out at the top of each page 
of output. 

A 0 allows the program to calculate N/A 

a shock wave while a 1 will cause 
the program to skip the shod calcu- 
lation. ICROSS is set equal to 
ICON(12) and ICROSS is used as a 
test for shock calculations. 
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Symbol 

ICON(13) 


ICON(14) 


ICON(15) 


ICON(16) 

IWALLi 


At nozzle 


Description 

Flag which tells program to use a vis- 
cous calculation in netting up the 
starting line* A zero does not read in 
viscous card and does not call subroutine 
VISCUS. A (l)reads in boundary layer 
information and calls subroutine VISCUS 
which calculates a boundary layer then 
merges the supersonic portion of the 
boundary layer into the starting line. 

ICON(14) is the number of Prandtl-Meyer 
rays to be used by the program in any 
expansions encountered in the solution. 

A zero will allow the program to calcu- 
late the number of rays. The program 
finds the limiting expansion angle then 
divides this angle into the number of 
rays to numerically integrate the 
Prandtl-Meyer expansion equation. 

ICON(15) is a two-digit number. The 
first digit is a flag which allows 
multiple cases to be run. When it is 
a 1 subroutine DRIVER reinitiates the 
program for the second case. A zero 
terminates the program after the first 
case has terminated. The second 
digit is used for reading in mesh con- 
trol criteria. If it is a 1, the mesh con- 
trol criteria are read in. If it is a 
zero, no mesh control is used in the 
solution. 

A flag which is set to a 1 will result in 
intermediate printout in the shock 
solution subroutines. There will be 
no intermediate printout if ICON(16) 
is a 0. This parameter is useful in 
program checkout. 

The type of boundary equation which is 
used in checking a characteristic 
point against problem boundaries. 

A 1 indicates a conic equation, a 2 
a polynomial equation and a 3 a free 
boundary equation. (See following page 
for description and example.) 


exit plane only. 
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A typical use of the conic equation is the throat of a nozzle. A sketch of 
a nozzle throat region shows the coefficients of circular throat. 



RC - radius of curvature of the 
circular arc of the throat 
RT - Throat radius 
XO - Axial distance from the origin 
of the coordinate system to the 
throat 

0 - Throat divergence angle correspond- 
ing to the maximum axial value for 
which the throat conic equation 
applies 


The conic equation for this case would have the following form 

A = -1 for an upper equation, +1 for a lower equation 
(-1 for this case) 

B = RC 2 - XO 2 
C = 2X0 
D = -1 

E = -(RC + RT) 

Xmax = RC sin 0 + XO 


An example of a free boundary is 



shown below. 

The free stream approach flow is 
inclined at 15 deg to the plume with 
a gamma (y) of 1.4, a Mach No. 
of 10., and a static pressure of 
0.1 PSFA. 


PINF 

E 

GAMMAINF 

MINF 

THETAINF 


0.1 PSFA 

0 (No pressure variation with axial distance) 
1.4 
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Symbol 

ITRANS 


WALLCO(5) 


XMAX 


RC 

RT 

THETA 


Description 


Indicates whether or not a discontinuity 
follows a particular boundary equation. 
ITRANS is used by the program in 
flaging calculations which will put a 
characteristic point at the exact loca- 
tion of the discontinuity (expansion or 
compression corner) so that a Prandtl- 
Meyer expansion or shock calculation 
may begin. 

The coefficients of the equations which 
describe the physical boundaries of the 
problems. These coefficients (A, B, 

C, D or E) and (PINF, GAMMAINF, 
MINF, THETINF or E) are used in 
either a conic or polynomial equation 
(A, B, C, D or E) or a free boundary 
equation (PINF, GAMMAINF, MINF, 
THETINF or E) depending on what was 
used for IWALL. 

XMAX is the maximum axial coordinate 
for which any one particular boundary 
equation is valid. XMAX is used by the 
program to test when the next boundary 
equation in the WALLCO(l 00, 6, 2) array 
should be used to check characteristic 
points against the physical boundaries 
of the problem. XMAX is stored in the 
6th position of the WALLCO array. 

RC, RT, THETA, and XO are used 
when it is desired to curve fit an upper 
or lower wall equation. RC is the 
radius of curvature of the nozzle throat 
and is used in calculating the equation 
of the nozzle throat. 

RT is the radius of the nozzle throat. 

RT is used in calculating the equation 
of the nozzle throat. 

The nozzle throat divergence angle. 
THETA is used to determine the maxi- 
mum axial coordinate for which the 
nozzle throat equation will be applied. 
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Units 

N/A 


The radial and 
axial components 
of the equations 
should be con- 
sistent with units 
for the particular 
problem being 
rum 


Same as above. 


Consistent with 
the units in 
which XIN and 
YIN are input 


Consistent with 
the units in 
which XIN and 
YIN are input 

Degrees 
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Symbol 


Description 


Units 


XO 


XIN 


YIN 

ALPHA 


UNITS 


If the nozzle throat is not at the origin 
then XO is the distance from the 
origin to the nozzle throat. 

The axial coordinate of a discrete 
point which is to be used in spline 
fitting the nozzle solid boundary. 


The radial coordinate of a wall point 
associated with a XIN coordinate. 

The gas name identifying the real gas 
properties on tape. This name is 
compared with the names of gas cases 
appearing on a master gas tape. When 
the program encounters the same case 
the gas data is read from the tape and 
stored for use by the program. If gas 
properties are to be read in from cards 
the only limit on the name is that it 
satisfy the format. This name is also 
written out at the top of the pages which 
contain the gas data written out. 

Either ENG or MKS, ENG can be used 
for cards only while MKS can be used 
for either cards or tape. If the gas 
properties are in MKS units, sub- 
routine GASRD changes the gas proper- 
ties to ENG units. 

ENG and MKS units are as follows: 


Consistent with 
the units in 
which XIN and 
YIN are input. 

The units of XIN 
and YIN should 
be the same as 
the units the pro- 
gram uses in the 
solution. 


N/A 


n/a 



ENG 

MKS 

Temperature 
Entropy 
Pressure 
Gas constant 

°R 

ft^/sec^°R 

psfa 

1545.5x32. 174 

°K 

cal/gram°K 

ATM 


MWT 

MWT (molecular weight) 


Note: These units are 
putting gas properties. 
ENG units internally. 

used only for in- 
Program uses 
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Symbol 


Description 


Units 


IOF 

IS 

OF RAT 

STAB 

IVTAB 

TAB(I, J, K,©) 

TAB (I, J, K,dJ 

TAB (I, J, K,©) 

TAB (I, J, K,@> 
TAB (I, J, K,©) 
CORLIP(2) 


The number of O/F ratio tables in the 
gas data. Used by the program for 
reading in gas properties from cards. 

The number of entropy cuts being con- 
sidered for each O/F ratio table. Also 
used for setting limits on calculations 
and reading in gas properties from cards. 

The O/F value for each O/F table. 

OF RAT is stored so that subroutine 
FABLE can locate the proper O/F tables 
to locate local gas properties in the flow 
field. 

An entropy value associated with a 
table of velocity cuts at an O/F 
ratio. Used in table lookup for local 
gas properties. 

The number of velocity cuts associated 
with a particular entropy cut. 

The mach no. associated with a particu- 
lar velocity cut and entropy for an O/F 
table. 

The gas constant associated with the 
mach number and entropy for an O/F 
table. 

GAMMA. The ratio of specific heats 
at a Mach no. and entropy for au O/F 
table. 

The static temperature at the Mach no. 
and entropy for an O/F table. 

The static pressure at the Mach no. and 
entropy for an O/F table 

The axial coordinate of the upper limit 
of the start line. CORLIP(2) is read in 

when it is desired for the program to 
set up the start line. CORLIP(2) is 
used to find the radial coordinate and 
flow angle at this axial coordinate. The 
resulting point is the first point on the 
start line. 


N/A 

n/a 


n/a 


ENG - — ^ — 
sec °R 

MKS 

gram K 

N/A 

n/a 

1545.4x32.174 

MWT 

MWT - MKS 
N/A 

°R - ENG 
°K - MKS 

PSFA - ENG 
ATM - MKS 

Consistant with 
boundary equa- 
tion. units. 
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Symbol 


Description 


Units 


CORLIP (6) 
CORLIP(4) 


CORLlP(5) 

CORLIP(8) 

CORLIP(9) 


STEP(6) 


PSI{? , M) 

PS 1(2, M) 

PS 1(3, M) 
PSI(4, M) 
PSI(5, M) 

PSI(6, M) 


The axial coordinate of the lower limit 
of the start line. Used for locating the 
lower point on the start line. 

The Mach no. or A/ A* for the start 
line location. If the starting line is set 
up using A/A*, subroutine AOASTR is 
called with A/A* and the Mach no. for 
the starting line is determined. The 
starting line Mach no. is used in sub- 
routine MASCON for determining a Mach 
no. distribution at the starting line. 

The entropy level of the starting line. 
Only used if the program sets up the 
starting line. 

The area of the nozzle throat. This is 
used for thrust coefficient calculations. 


Consistent with 
boundary equa- 
tion units. 

N/A 



Consistant with 
boundary equa- 
tions. 


A constant O/F ratio across the starting N/A 

line. This value is used only if vhe pro- 
gram is setting up the starting line. 

Note: The CORLIP values are used only 
if the program is to set up the starting 
line. 

The minimum percent change in pitot total N/A 
pressure across a shock for discontinuing 
shock calculation. This is used to test 
shock strength. If the percent change from 
oblique shock theory is less than Step (6), 
then no iteration is performed by subroutine 
SHOCK, thereby keeping the shock strength 
the same. If nothing is read into the pro- 
gram, Step (6) is set to .001. 


For an input starr line this is the radial 
coordinate for a particular starting 
line point. 


Consistent with 
boundary equa- 
tion. 


The axial coordinate of each start line 
point. 

The Mach no. at each start line point. 

The flow angle at each start line point. 

The entropy level at each starting line 
point. 

The shock angle downstream of the 
shock point on the start line. This 
value is input on first card only for a 
reflected shock case. 


Consistent with 
boundary equa- 
tions. 

N/A 

Degrees 

ft 2 /sec 2o R 

Degrees 


3-21 


LOCKHEED - HUNTSVtLIE RESEARCH & ENGINEERING CEN . LR 



LMSC/HREC D 162220 -IV -A 


Symbol 


Description 


Units 


PSI(8, M) 
CUTDAT(l) 

CUTDAT(2) 


CUTDAT(3) 


CUT DAT (4) 


CUTDAT(5) 

CUTDAT(6) 


The value of the O/F ratio at each 
starting line point* 

The radial coordinate defining upper 
cutoff of the characteristic solution* 


The axial coordinate defining upper 
cutoff* If a characteristic point is 
calculated whose axial coordinate 
falls behind this coordinate the char- 
acteristic line is terminated and the 
next line started. 

The angle the upper cutoff line makes 
with the horizontal. This line initiates 
from the point formed by CUTDAT(l) 
and CUTDAT(2)* If a characteristic 
point is calculated which falls above 
this line the characteristic line is ter- 
minated and a new line begun. 

The radial coordinate defining down- 
stream cutoff. If a characteristic 
point is calculated whose radial 
coordinate falls below this coordinate 
the line is terminated and a new char- 
acteristic line started. 

The axial coordinate defining down- 
stream. cutoff. 


The angle the downstream cutoff line 
makes with the horizontal. The cut- 
off line initiates from the point formed 
by CUTDAT(4) and CUTDAT(5). If a 
characteristic point is calculated which 
falls beyond (and direct) this line the 
characteristic line is terminated and a 
new line begun. 

Note: The problem is terminated when 
a characteristic line is started beyond 
the limits set by CUTDAT 1 s 1-6. 


N/A 


Consistent with 
boundary equa- 
tions . 

Consistent with 
boundary equa- 
tions. 


Degrees 


Consistent with 
boundary equa- 
tions 


Consistent with 
boundary equa- 
tions. 

Degrees 
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Example of Cutoff Limits 




© Lines Terminated due to Cutoff Limits 
and Free (Pressure) Boundary 


- CLTDAT(i) 

X l - CUTDAT(2) 

Q l - CUTDAT(3) 

R 2 - CUTDAT (4) 

X 2 - CUTDAT(&) 

e 2 - CUTDAT(6) 

0 - Line Terminated 
because of Cutoff 
Limits 

0 - Line Terminated 
because of Free 
Boundary 

• Characteristic 

Data Points along 
left -running char- 
acteristic lines 
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Symbol 

DIV(l) 

DIV(2) 


DIV(3) 

DIV(4) 

DIV(5) 

DIV(6) 


Description 

The maximum change in absolute length 
(axial or radial) between either of the 
two base points and a new characteris- 
tic point. If the change in distance is 
greater than this maximum a new point 
is added by the program. 

The minimum change in absolute length 
(axial or radial) between either of the 
two base points and a new characteris- 
tic point. If any change in length is 
less than this minimum the new char- 
acteristic point will be deleted. 


The maximum change in Mach no. be- 
tween either of the two base points and 
a newly calculated base point. If the 
change in Mach no. exceeds this maxi- 
mum a new point will be inserted which 
will satisfy the criteria. 

The minimum change in Mach no. be- 
tween either of the two base points and 
a newly calculated base point. If the 
change in Mach no. is not greater *han 
this value the point will be deleted. 


The maximum change in flow angle be- 
tween either of the two base points and 
the newly calculated point. If the change 
in flow angle is greater than the maximum 
a new point is added which will satisfy the 
criteria. 

The minimum change in flow angle be- 
tween either of the two base points and 
the new characteristic point. If the 
change in flow angle is less than this 
minimum then the point is deleted. 


Note: Mesh control explained in detail 
in Appendix. 
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Consistent with 
boundary equa- 
tions. 


Consistent with 
boundary equa- 
tions. 


N/A 
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Symbol 


Description 


Units 


NPOWER 


XL 


NBLTS 

CU 


The reciprocal of the exponent of the N/A 

velocity profile in the boundary layer. 

This exponent is used only if the Reynolds 
number calculated for the case at the exit 
conditions at the wall indicates turbulent flow. 

If the flow is found to be laminar NPOWER 
is set to 2 by the program. 


A characteristic length which will be 
used to calculate the Reynolds no. and 
the boundary layer thickness. This 
value is usually taken as the nozzle wall 
length. 

The number of starting line points which 
are desired in the boundary layer. 

A conversion factor for mixed units of 
length. This conversion factor should 
be such that XL(ft) will be converted 
into the same units used in the boundary 
equation. If inches are used in boundary 
equations, then CU will be 12. If boundary 
equations are in ft., CU is 1. CU is used 
in calculating the boundary layer thickness. 


Ft 


N/A 

inches /ft or 
as required 
by boundary 
equation units. 


NOTE: A detailed description of the boundary layer option 
is explained in .the Appendix. 


STEP(3) 


TINF 


Point insert criteria. Step(3) is the maxi- consistent 
mum change in axial distance between two with flow- 
consecutive characteristic lines on a lower field units 
wall. If this value is exceeded then a new 
left running characteristic line is started at 
the lower wall. 


Free stream static temperature. This °R 

temperature is used for the oblique shock 
calculation of the plume free boundary 
when IWALL on card 3 is set to a 6. 
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3.2 PROBLEMS COMMONLY ENCOUNTERED 


This section is intended to aid the user in utilizing the program and 
avoiding some common problems. The comments on mesh control and start- 
ing line points are only initial guidelines. A "feeling" for the starting line 
points and mesh control criteria will result as the user becomes familiar 
with the operation of the program. 


The following is a list of hints to the user: 

• If an automatic restart is not desired when a shock intersects the 
axis, set ICON(6) to zero in Cols 28,29 and 30 on Card 2. (Refer 
to Section 3.1.2 and the Appendix.) 

• If it is desired to keep the type of printed output the same through- 
out the entire plume use only Cols. 39 and 40 of Card 2 (ICON(8)). 

• For inputting a start line for a highly expanded plume the start 
line points should be more highly concentrated near the nozzle wall. 

• When setting up a nozzle wall with discrete points the first spline 
fit point must not be a part of the nozzle throat equation (i.e., the 
first input point should be beyond the maximum axial coordinate 
for which the throat equation applies). 

• In general, the mesh control criteria can be relaxed as the 
number of start line points increases. 

• For highly expanded nozzles where a very large plume is to be 
generated an approximate mesh control criteria for the maximum 
Mach quadrilateral size is 20 exit diameters. Corresponding 
values for Mach number and flow angle criteria would be about 

10 and 20° respectively. 

• If problems in the solution are encountered due to characteristic 
lines diverging, decrease the maximum mesh control criteria so 
that more lines may be added. 

• When the messages "negative velocities encountered" or "subsonic 
Mach number encountered," first check the input data. If no errors 
are noted in the input data decrease the mesh control criteria, or 
alter the number of start line points. 
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3.3 BASIC LIST OF SUBROUTINES, SIMPLE PROGRAM FLOWCHART 
AND COMMON BLOCK VARIABLES 

• General Flowchart 

• Lists of Subroutines according to function within the program 

• Definition of all common block variables 
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Fig. 3-1 - General Flow Chart 
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3.3.1 Program Routines 

The following subsection contains a list of the routines according to 
their function performed within the program. Included are: 

• General Flow Properties Routines 

• Characteristic Solution Routines 

• Shock Routines 

• Logic Control 

• Tape Manipulation Routines 

• Start Line Setup Routines 

• Problem Limit and Boundary Routines 

• Input/Output Routines 

• Utility Routines 
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GENERAL FLOW PROPERTIES 
ROUTINES 


Routine 


Function 


EMOFP 

EMOFV 

FABLE 

FNEWTN 

HYPER 

POFEM 

PRANDT 


RGMOFP 

RGVOFM 

RHOFEM 

TABLE 

THETPM 

TOFEM 

TOFV 

UOFEM 

UOFV 

VOFEM 


finds Mach number as a function of local static pressure 
and entropy. 

finds Mach number as a function of velocity. 

finds local gas properties as a function of O/F ratio 
with the aid of (TABLE) 

finds Newtonian impact pressure at the plume boundary. 

calculates the balanced pressure at the intersection of a 
solid boundary and pressure boundary. 

finds static pressure as a function of Mach number and 
entropy. 

finds the Prandtl-Meyer expansion angle for a given 
boundary angle and divides the angle into a series of 
expansion "rays". Then the flow properties are deter- 
mined for each ’'ray". 

finds Mach number as a function of static pressure, 
entropy and O/F ratios. 

finds velocity as a function of Mach number, entropy 
and O/F ratio. 

finds density as a function of Mach number and entropy. 

finds local gas properties for an entropy and velocity 
within an O/F table. Uses tables of gas properties in- 
put to the program. 

performs a numerical integration to calculate properties 
through a Prandtl-Meyer expans ion. 

finds static temperature as a function of Mach number, 
finds static temperature as a function of velocity, 
finds Mach angle as a function of Mach number, 
finds Mach angle as a function of velocity, 
finds velocity as a function of Mach number. 
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Routine 


MOCSOL 

MONO 

ROTERM 


CHARACTERISTIC SOLUTION ROUTINES 

Function 


provides all two-dimensional and axisymmetric methods 
of characteristic solutions. 

determines if the characteristic point locations are mono- 
tonic along left or right running lines. 

finds the geometrical factor used in the axisymmetric 
term of the compatability equation and as an interpola- 
tion parameter. 
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SHOCK ROUTINES 


Routine 


Routine 


CONIC 

DELTAF 

ENTROP 

ESHOCK 

OVEREX 

SHOCK 

TURN 

WEAK 


solves for flow properties immediately downstream of 
a shock initiating from conic surface immersed in the 
flow. 

computes the turning angle through an oblique shock 
knowing the shock angle and the upstream Mach no. 

finds entropy rise across a shock as a function of shock 
angle and upstream Mach number. 

uses an iterative solution to perform equilibrium shock 
calculations. 

finds the shock angle at the nozzle lip for over expanded 
flow. 

iteratively adjusts the shock strength in order to satisfy 
the oblique shock relations and the flow field properties 
s imultaneously. 

solves for a shock wave which has a known turning 
angle. 

finds entropy and velocity downstream of the shock. 
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Routine 

DRIVER 

KIKOFF 

PHASE1 

BLCKBX 


LOGIC CONTROL 

Function 


provides highest order control for program execution. 
The initialization and logic subroutines are called 
from here. 

allows control to return to DRIVER if an error is en- 
countered in the calculations, 

provides the necessary logic to employ the proper cal- 
culations at the proper time in order to describe the 
entire characteristic mesh, 

provides controlling logic for automatic restart of a 
reflected shock case. 
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TAPE MANIPULATIONS ROUTINES 


Routine 


Function 


GASRD 

GASTAP 

IDTAPE 

JWBTSS 

OUT BIN 
READB 
READF 
TAPMOV 


reads gas properties from cards and converts gas 
properties from MKS to English (ENG) units if neces- 
sary. 

reads real gas properties from tape generated by the 
NASA Lewis Thermochemical Data program and writes 
data on output tape. 

writes gas properties read from cards on flow field 
tape. 

finds flow properties at es ;h point along the starting 
line for the automatic restart of the reflected shock 
case. 

writes the calculated characteristic point data on the 
binary output tape. 

reads the position of the boundary points from the flow- 
field tape. 

reads one characteristic line from the flowfield tape and 
saves flow properties at each point on the line. 

moves the flowfield tape through the gas properties to 
the beginning of the flowfield information. 


3-34 


LOCKHEED -HUNTSVILLE RESEARCH & ENGINEERING CENTER 



LMSC/HREC D162220 -IV 


START LINE SETUP ROUTINES 


Routine 


Function 


AOASTR 

LIPIN 

MASCON 

SETUP 

VISCUS 

WOFA 


finds the Mach no. corresponding to a given area ratio. 

calculates information for starting line points setup 
by the program. 

calculates Mach no. distribution at any area downstream 
of the nozzle throat while conserving total mass flow. 

locates and distributes start line points for the automatic 
restart of the reflected shock case. 

calculates the boundary layer thickness at the nozzle 
exit and adjusts the starting line to consider the boundary 
layer. 

finds weight flow per unit area as a function of Mach 
number. 
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PROBLEM LIMITS & BOUNDARY ROUTINES 


Routine 

Function 

BOUND 

finds the radial coordinate and flow angle for a given 
axial coordinate on an upper or lower solid boundary. 

ITERM 

tests each characteristic point to determine if it lies 
within problem limits. 

LIMITS 

determines which boundary equation to use in checking 
on problem limits. 

SOLVE 

finds the coefficients of the nozzle wall cubic equation 
described by points input to the program. 

SWITCH 

transfers upper and lower boundary equations for 
automatic reflected shock restart case. 

THROAT 

computes the nozzle throat equation for a nozzle wall 
described by discrete points. 
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Routine 

ERRORS 

OUT 

PAGE 

PLUMIN 

PLMOUT 


INPUT/OUTPUT ROUTINES 

Function 


writes out various messages for errors which may 
occur during program execution. 

writes the calculated data at characteristic points along 
with the title and heading. 

page ejects and writes header comments and page num- 
ber on each page of printout. 

reads in the input data necessary to perform the method 
of characteristics solution. 

prints data read by (PLUMIN). 
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MISCELLANEOUS ROUTINES 


Routine 

Function 

ICHECK 

this routine determines whether a characteristic point 
is added or deleted depending on v he mesh control 
criteria. 

INITP 

initializes values for various control parameters. 

INRSCT 

find intersection of two straight lines. 

ITSUB 

controls the iteration for any set of equations that are 
expressed as a function of one variable. 

MASSCK 

keeps a running check on the mass flow under each 
point along a characteristic line. 

PRANDT 

finds the Prandtl-Meyer expansion angle for a given 
boundary angle and divides. 

SAVER 

saves data at nozzle lip. 

THRUST 

computes the vacuum thrust produced by a two- 
dimensional or axisymmetric nozzle. 
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3.3.2 Common Variables 


This subsection contains a list of the common block variables along with 
a description of each variable and the routines using each common block. 


Common 

Name Variables Definition Routines Using 


/CARDT/ icard 


/CONTRL/ IRUN(IO) 
ICON(16) 


If ICARD is a 1 subroutine 
IDT APE will write gas data 
on output tape. If ICARD is 
not a one subroutine GASTAP 
will read the master gas tape 
or write the gas data on the 
output tape. (ICARD is used 
in automatic shock reflection. ) 

Control flags used by the pro- 
gram coring the plume genera- 
tion. 

Input controls used by the 
program. 


/CRITER/ CONVRG(IO) Array of convergence criteria. 

ITLIM(IO) Array of iteration limits. 


BLCKBX 

GASRD 

DRIVER 


DRIVER 

ESHOCK 

FABLE 

GASRD 

GASTAP 

IDTAPE 

INITP 

KIKOFF 

LIMITS 

LIPIN 

MASCON 

MASSCK 

MOCSOL 

OUT 

PAGE 

PHASE1 

PLMOUT 

PLUMIN 

PRANDT 

RGMOFP 

RGVOFM 

SHOCK 

TABLE 

THRUST 

TURN 

VISCUS 

DRIVER 

INITP 

MOCSOL 

PRANDT 

SHOCK 

TURN 
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Common 

Name Variables Definition Routines Using 


/CUTFO/ CUTDAT(6) cutoff limitations (problem 

limits) 


/DATAR/ PHO(8, 100, 2) 

PHO(100, 2) 

APHO(100, 2) 

OFRAT(IO) 
STAB(10, 2) 

IVTAB(10, 2) 

TAB(10,2,13; 

5) 

WALLCO(100, 

6 , 2 ) 

IWALL(100,2) 

ITRANS(100, 

2 ) 

IEQNOW (2) 


Local gas properties array 
for all characteristic points 
on the line the program is 
calculating and the previous 
characteristic line. 

Type of characteristic point 
for each point on the new 
characteristic line and the 
previous one. 

Array which stores the Mach 
no. for et.ch point on the new 
characteristic line and the 
previous one. 

O/F ratio for each O/F table. 

Entropy values for each O/F 
ratio table. 

Number of velocity cuts in 
each entropy table for each 
O/F ratio. 

Local gas properties for each 
velocity cut in each entropy 
table for each O/F ratio. 

Wall coefficients and maxi- 
mum X for each upper and 
lower wall equation. 

Type of each upper and lower 
wall equation. 

Contains flag for each wall 
equation telling program 
whether an expansion corner, 
compression corner or no 
discontinuity follows the 
equation. 

Counter on the upper and 
lower boundary equations. 


DRIVER 

ITERM 

PLMOJT 

FLUMIN 

BOUND 

DRIVER 

FABLE 

FNEWTN 

GASRD 

GASTAP 

HYPER 

ICHECK 

IDTAP 

INITP 

ITERM 

LIMITS 

MASSCK 

MOCSOL 

MONO 

OUT 

OUT BIN 

OVEREX 

PHASE 1 

PLMOUT 

PLUMIN 

PRANDT 

RGMOFP 

RGVOFM 

SHOCK 

SOLVE 

TABLE 

THROAT 

THRUST 
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Common 

Name Variables Definition 


Routines Using 


/DEL/ IDEL 

DECRE 

/DIVCRl/ DIV(6) 

/fab/ t 

p 

SP 

/FLAG/ iflag 

JFLAG 

KFLAG 

/FORCE/ FORCEX 

FORCEY 

TORQZ 

EMASS 


Flag used for printing out 
which of the mesh control 
criteria a point has not satis- 
fied. 

The distance, Mach no. or 
flow angle for the point which 
did not satisfy the mesh con- 
trol criteria. 

Mesh control parameters 
which control mesh size by 
inserting or deleting points. 


Local temperature at point 
of interest. 

Local pressure at point of 
interest 

Local entropy at point of 
interest 

Determines if gas properties 
have been read into program; 

1 - not read in - 2 - have been 
read in. 

Indicates if last point was in 
plume; 2 - last point in plume; 
1 - last point not in plume. 

Indicates if present point is 
in plume; 0 - not in plume; 

1 - is in plume. 

Total axial force on nozzle 
up to the last calculated 
characteristic line. 

Total radial force on nozzle 
up to the last calculated 
characteristic line. 

Total moment about the origin 
due to FORCEX and FORCEY. 

Total mass flow passing 
through the nozzle. 


DRIVER 

ICHECK 

PLMOUT 

PLUMIN 

FABLE 

TABLE 


BLCKBX 

JWBTSS 


DRIVER 

MASSCK 

THRUST 
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Common 


Name 

Variables 

Definition 

Routines Using 

/GASCON/ 

R 

Local gas constant in plume 

DELTAF 


GAMMA 

Local isentrepic exponent in 
plume. 

DRIVER 

EMOFP 

EMOFV 


TO 

Total temperature. 

ENTROP 

ESHOCK 


PO 

Total pressure. 

FABLE 


ASTER 

Flag which indicates whether 
gas properties are in table or 
not. 

GASRD 

MOCSOL 

OUT 


PHASE1 

PLMOUT 

PLUMIN 

POFEM 

PRANDT 

RGMOFP 

RHOFEM 

SHOCK 

TABLE 

THETPM 

TOFEM 

TOFV 

VISCUS 

VOFEM 

WEAK 

WOFA 

/HEAD/ HEADER(36) Array containing the title of DRIVER 

the MOC run- GASTAP 

IDTAP 

INITP 

OUT 

PAGE 

PI MOUT 

PLUMIN 


/INPUT/ PSI(8, 100) Local gas properties array for DRIVER 

each point on the starting line. LIPIN 

MASSCK 

PH AS El 

PLMOUT 

PLUMIN 

THRUST 

VISCUS 
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Name 

/IRECD/ 

/IRFEC/ 

/ITAPE/ 

/points/ 


/QUIT/ 

/SHCKPT/ 


P 

LMSC/HREC D1 62220 -IV 


Variables 


Definition 


Routines Using 


IREC A counter on the number of 

files written on the binary- 
output tape. 

IFILE Flag which tells the program 

where to begin writing on the 
binary output tape. 

IGORN Flag which causes a message 

to be written out when the 
shock intersects the axis. 


ISTARO(50, 3) Contains the point number of 
the beginning of a new char- 
acteristic line which has been 
added due to mesh control. 
Also contains the point num- 
bers of the two points used to 
generate the new point. 

IREST The line number identifying a 

characteristic which has been 
added due to mesh control. 


XIN(100, 2) Axial location of points input 
to describe the upoer or 
lower nozzle wall points. 

YIN(100, 2) Radial location of points input 
to describe the upper or lower 
nozzle wall points. 


RC 


Radius of curvature of the 
nozzle throat. 


RT 


Radius of nozzle throat. 


DRIVER 

TAPMOV 


DRIVER 

SHOCK 


PH AS El 
OUTBIN 


PLUMIN 

SOLVE 

THROAT 


THETA 

XO 

ISTOP 

PMX(8) 


Nozzle throat divergence angle. 
Axial coordinate shift. 


0 - one case is to be run. 

1 - more than one case is 

to be run. 

Array which contains flow 
properties at the nozzle lip 
point or last calculated 
nozzle Weill point (used for 
automatic shock reflections). 


DRIVER 

PLUMIN 


BLCKBX 

SAVER 

DRIVER 
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Common 

Name 

/STEPC/ 

/TAPEFO/ 

/UNITSS/ 

/WHICH/ 


/WTSAVF/ 

/XSICOM/ 


Variables 

Definition 

Routines Using 

STEP(IO) 

Array containing floating 
point information used as 
controls by the program. 

DRIVER 

HYPER 

INITP 

MASSCK 

PHASE 1 

PLUMIN 

PRANDT 

THETPM 

THRUST 

ALPHA(4) 

Identification of propellant 
case. 

DRIVER 

FABLE 

IOF 

Number of O/F ratio tables. 

GASRD 

GASTAP 

IS 

Number of entropy cuts for 
each O/F ratio table. 

IDTAP 

PLMOUT 

RGMOFP 

RGVOFM 

UNITS 

The type of units the gas 
properties were read in as. 

DRIVER 

GASRD 

ENG 

Gas properties are in English 
units. 


EMKS 

Gas properties are in MKS 
units. 


KAY 

The line on which the shock 
wave has intersected the 
axis. 

BLCKBX 
DRIVER 
PLUMiN 
PHASE 1 

IWHARE 

The boundary equation at 
which shock restart calcula- 
tions are to begin. 


ILINE 

The type of printed output de- 
sired when the program begins 
using a free boundary equation 
for problem units. 


WTSAV(IOO) 

Total mass flow passing under 
each point on a characteristic 
line. 

MASSCK 

OUTBIN 

XSI(10, 2, 
13,2) 

Weighing factors for interpo- 
lating between velocity cuts 
over pressure and temperature. 
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Common 

Name Variables Definition Routines Using 


/FINITE/ IFREE 


TINF 


A flag which is used to deter- ESHOCK 
mine whether a free boundary- FABLE 
is to be calculated using oblique FNEWTON 
shock theory. PHASE 1 

TURN 


Freestream static temperature 


OF 


Not presently used 
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3.4 DESCRIPTION AND FLOW CHARTS OF INDIVIDUAL SUBROUTINES 

FUNCTION NAME: AOASTR 

DESCRIPTION 

This function finds the Mach number corresponding to a given area 
ratio by one -dimensional theory. Real gas effects are considered in this 
calculation. 

CALLING SEQUENCE 

EM = AOASTR (OF, S, AOA) 

where (EM) is the Mach number which exists, one dimensionally, at an area 
ratio of (AOA), an entropy (S), and at an O/F ratio (OF). 

UTILITY ROUTINES AND COMMON REFERENCES 

COMMON- None 

ERRORS 

FABLE 

ITSUB 

RGVOFM 

WOFA 

METHOD OF SOLUTION 

The weight flow per unit area at Mach one is evaluated. An initial 
guess for the desired Mach number is made and ITSUB is initialized. An 
iterative solution of the equation FOFEM = AOA - WOFAl/WOFA(EM), 
driving FOFEM to zero, is performed with the aid of ITSUB. 
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SUBROUTINE NAME: BLCKBX 


DESCRIPTION 

This routine contains the controlling logic for automatically restarting 
the method of characteristics solution for a reflected shock case. 

CALLING SEQUENCE 
CaU BLCKBX 

UTILITY ROUTINES AND COMMON REFERENCES 


COMMON/CARDT/ 

plmout 

common/divcri/ 

idtape 

common/cutfo/ 

gastajp 

COMMON/ST EPC/ 

TAPMOV 

common/input/ 

FABLE 

common/flag/ 

EMOFV 

common/shckpt/ 

SETUP 

common/dat ar/ 

JWBTSS 

common/cont rl/ 

SWITCH 

common/which/ 



METHOD OF SOLUTION 

Not applicable. 
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SUBROUTINE NAME: BOUND 
DESCRIPTION 

This subroutine £inds the radial coordinate and flow angle (radians) for 
a given axial coordinate on an upper or lower solid boundary. 

CALLING SEQUENCE 

CALL BOUND (R, X, THETA, ITYPE) 

where (R) is the radial coordinate, (X) is the known axial coordinate, THETA 
is the wall angle and ITYPE indicates whether upper or lower boundary equa- 
tions are to be used. 

UTILITY ROUTINES AND COMMON REFERENCES 

COMMON/D AT AR / 

UTILITY - None 

METHOD OF SOLUTION 

The block common region DATAR contains boundary equations neces- 
sary to evaluate R, THETA. Two types of equations are used; 

r = a [ *Jb + cx + dx 2 + e ] CONIC TYPE 1 

and 

r = ax 4 + bx 3 + cx 2 + dx + e POLYNOMLVL TYPE 2 

The input fixed point variable ITYPE has a one or a two in the units position 
which selects the upper (2) or lower (1) coefficients and control information. 
IEQNOW contains the number of the equation to be used. 
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SUBROUTINE NAME: CONE 
DESCRIPTION 

This subroutine solves for the shockwave downstream properties formed 
at a right circular cone immersed axisymmetric supersonic flow. 

CALLING SEQUENCE 

CaU CONE (OF, SU, QU, ALPHA, SD, QD, EPS) 

where the input properties are (OF, SU, QU) the upstream O/F ratio, entropy, 
and velocity and (ALPHA) is the turning angle. The subroutine returns with 
(SD, QD) the downstream entropy and velocity and (EPS) the shock angle. 

UTILITY ROUTINES AND COMMON REFERENCES 

COMMON - None 

FABLE 

TOFV 

ESHOCK 

EMOFV 

ITSUB 

ERRORS 

METHOD OF SOLUTION 

The oblique shock relations and characteristic relations are solved in 
an iterative manner to obtain the downstream flow properties. 


3-49 


LOCKHEED • HUNTSVILLE RESEARCH & ENGINEERING CENTER 



LMSC/HREC D162220-IV 


FUNCTION NAME: DELTAF 
DESCRIPTION 

This function computes the turning angle through an oblique shock wave 
knowing the shock angle and the upstream Mach number. 

CALLING SEQUENCE 

DELTA = DELTAF (EPS, EM) 

where (DELTA) the turning angle is found from the shock angle (EPS) and 
the upstream Mach number (EM). NOTE: The appropriate values of gas 
properties must be stored in common upon entry to this routine. 

UTILITY ROUTINES AND COMMON REFERENCES 

COMMON/C ASCON/ 

UTILITY - None 

METHOD OF SOLUTION 

The oblique shock relations are solved for the turning angle using the 
relations; 

6 = e • tan ' 1 + M-t)| 
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SUBROUTINE NAME: DRIVER 
DESCRIPTION 

Driver provides the highest order control for program execution. 
Initialization and logic subroutines are called from here. Most all the 
common storage needed in the remainder of the program is specified here. 

CALLING SEQUENCE 

CALL DRIVER (K) 

where (K) is a control constant indicating whether or not errors exist in the 
execution of the program. (K= 1 for a detected error, K = 0 for no errors.) 


UTILITY ROUTINES AND COMMON REFERENCES 


common/control/ 

COMMON/QUIT/ 

common/criter/ 

COMMON/XSICOM/ 

common/cutfo/ 

COMMON/CARDT/ 

common/ d at ar/ 

common/flag/ 

common/divcri/ 

common/shckpt/ 

common/force/ 

common/irecd/ 

COMMON/GASCON/ 

common/irfec/ 

common/head/ 

common/unitss/ 

common/input/ 

INITP 

common/stepc/ 

PLUMIN 

common/tapefo/ 

PHASE 1 


BLCKBX 

METHOD OF SOLUTION 


Not applicable for this subroutine. 
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FUNCTION NAME : EMOFP 
DESCRIPTION 

This routine computes the local Mach number as a function of local 
pressure (static) and the local entropy value. 

CALLING SEQUENCE 


EM = EMOFP (P.S) 

where (EM) is the resultant Mach number found from the pressure (P) and 
entropy (S). NOTE : The appropriate values of the gas properties must be 
stored in common upon entry to this routine. 

UTILITY ROUTINES AND COMMON REFERENCES 

C OMM ON/GASCON/ 

UTILITY - None 

METHOD OF SOLUTION 


Perfect gas relationships (thermally perfect) are used to find the 
Mach number. 
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FUNCTION NAME : EMOFV 
DESCRIPTION 

This routine finds Mach number as a function of velocity. 

CALLING SEQUFNCE 

EM 9 EMOFV(V) 

where (EM) is the local Mach number which is found as a function of (V) the 
local velocity. NOTE : The appropriate values of the gas properties must 
be stored in common upon entry to this routine. 

UTILITY ROUTINES AND COMMON REFERENCES 

common/gascon/ 

TOFV 

METHOD OF SOLUTION 

Thermally perfect gas relationships are used to find the Mach number. 

m -V(t • 
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FUNCTION NAME ; ENTROP 
DESCRIPTION 

This routine utilizes the oblique shock relations to find the entropy- 
rise across a shock as a function of the shock angle and the upstream Mach 
number. 

CALLING SEQUENCE 

SD = ENTROP (EPS, EMU) 

where (SD) is the entropy rise across the shock and is a function of the 
shock angle (EPS) and the upstream Mach number (EMU). NOTE : The 
appropriate values of the gas properties must be stored in common upon 
entry to this routine. 

UTILITY ROUTINES AN ..OMMON REFERENCES 

common/gascon/ 

UTILITY - None 
METHOD OF SOLUTION 

The oblique shock relations are employed to find the entropy rise 
across the shock. 

* . ^ - ( Z. r D»] + r fc fa^]} 


3-54 


LOCKHEED -HUNTSVILLE RESEARCH & ENGINEERING CENTER 




p 


p 


LMSC/HREC D 162220 -dV 


SUBROUTINE NAME : ERRORS 
DESCRIPTION 

This subroutine contains print messages for various errors which may 
occur. This is an open ended routine in that it can easily be extended to 
handle more print messages. 

CALLING SEQUENCE 

CALL ERRORS (I) 

where (I) selects the message to be printed. 

UTILITY ROUTINES AND COMMON REFERENCES 
None 

METHOD OF SOLUTION 
Not applicable. 
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SUBROUTINE NAME: ESHOCK 
DESCRIPTION 

This subroutine employs an iterative solution to perform the equili- 
brium shock calculations for a real or ideal gas. The real and ideal gas 
calculations are similar, the difference being that an ideal gas case con- 
verges on the first iteration. 

CALLING SEQUENCE 

CALL ESHOCK (OF, SI, VI, EP, DELTA, S2, V2) 

where the input properties are (OF) the upstream O/F ratio, (SI, VI) thft 
upstream entropy and \ elocity and (EP) the shock angle. The subroutine 
returns with (DELTA), the tur ning angle and (S2, V2), the downstream 
entropy and velocity. 

UTILITY ROUTINES AND COMMON REFERENCES 

COMMON/CONTRL/ POFEM 

COMMON/GASCON/ DELTAF 

COMMON/FINITE/ ENTROP 

EMOFV RHOFEM 

FABLE WEAK 

METHOD OF SOLUTION 

The continuity equation coupled with the equations for conservation of 
normal and tangential momentum are solved in an iterative manner utilizing 
thermochemical property data to satisfy the conservation of energy. This 
set of 4 equations is expressed in terms of the 4 unknown quantities: 

€ - shock angle 

6 - turning angle 

S^- entropy downstream of shock 

Vg" velocity downstream of shock 
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SUBROUTINE NAME: FABLE 
DESCRIPTION 

This subroutine utilizes real or ideal gas information obtained from 
the flowfield tape and a local O/F ratio to call subroutine TABLE to calcu- 
late thermodynamic gas properties locally in the flow. 

CALLING SEQUENCE 

CALL FABLE (OF.SS.V) 

where (OF) is the local O/F ratio, (SS) is the local entropy and (V) is the 
local velocity. 

UTILITY ROUTINES AND COMMON BLOCKS 

C OMM ON/GAS C ON/ 

COMMON/DATAR/ 

COMMON/FAB/ 

COMMON/CONTROL/ 

COMMON/FINITE/ 

TABLE 

METHOD OF SOLUTION 

The routine is entered with the local O/F ratio (OF), entropy (SS), 
and velocity (V). The local ratio is used to determine which set of thermo- 
dynamic tables that subroutine TABLE should use to perform table lookup 
of the local thermodynamic gas properties. Subroutine FABLE then uses 
the local thermodynamic gas properties obtained from TABLE to perform 
an interpolation between the O/F tables based on the local O/F ratio. 
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FUNCTION NAME : FNEWTN 
DESCRIPTION 

This function solves for the Newtonian impact pressure along the plume 
boundary. The calculation is applicable for hypersonic free stream velocities 
or quiescent conditions (i.e., = 0). There is also an option whereby im- 

pact pressure at a free boundary is solved using oblique shock theory. 

CALLING SEQUENCE 

P„, = FNEWTN (THETA3, X, ITYPE1) 

1M 

where is the hypersonic Newtonian impact pressure at the plume 

boundary, (THETA3) is the local flow angle at the boundary, (X) is the 
axial coordinate of the boundary point, and (ITYPE1) designates upper of 
lower boundary. (1 - lower, 2 - upper) 

UTILITY ROUTINES AND COMMON REFERENCES 
COMMON/GASCON/ 

C OMM ON/C ON TR OL/ 

COMMON/FINITE 

TURN POFEM VOFEM 

JOFEM EMOFV 

METHOD OF SOLUTION 

The block common region - WALLCO contains the necessary informa- 
tion to evaluate the free stream gas properties at the plume boundary point. 
The impact pressure is then calculated using the following equation 

P = P (1 + eX) fl + y M 2 sin 2 (0„ - 9 H 
oo L oo oo B oo J 

If the oblique shock option is used the free stream gas properties and 
impact pressure are calculated using subroutine TURN. 
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SUBROUTINE NAME: GASRD 
DESCRIPTION; 

This subroutine reads in the gas properties. These properties may be 
real or ideal and read in via cards or tape. The routine also converts input 
gas properties from MKS units to English (ENG) units if necessary. 

CALLING SEQUENCE 

CALL GASRD 

UTILITY ROUTINES AND COMMON REFERENCES 

common/xsicom/ 

common/contrl/ 

common/datar/ 

COMMON/GASCON/ 

COMMON/TAPEFO/ 

COMMON/UNITSS/ 

COMMON/CARDT/ 

GAS TAP 
IDTAPE 

METHOD OF SOLUTION 

The gas name (ALPHA(I)), type units, number of O/F tables, and num- 
ber of entropy cuts are read in from an input card. If the gas properties are 
on cards, this subroutine reads the cards. If the gas properties are on tape, 
control of the reading of properties is given to GASTAP. In either case, t’.e 
properties are converted from MKS to English (ENG) units by this subroutine 
if necessary. 
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SUBROUTINE NAME: GASTAP 
DESCRIPTION 

GASTAP retrieves the real gas properties generated by the NASA 
LEWIS THERMOCHEMICAL DATA program and provides instructions for 
writing this data on the MOC flow field tape. 

CALLING SEQUENCE 

CALL GASTAP 

UTILITY ROUTINES AND COMMON REFERENCES 

COMMON/CONTRL/ 

COMMON/DATAR/ 

common/head/ 

common/tapefo/ 

ERRORS 

METHOD OF SOLUTION 

The gas name, (ALPHA(I)), specified on the input data is compared 
with available cases on the thermochemical data tape until a match is found. 
This particular case is then read, stored core in and written on the MOC 
flow field tape. 
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SUBROUTINE NAME: HYPER 
DESCRIPTION 

This subroutine calculates the balanced pressure at a corner point 
(i. e. , at the intersection of a solid boundary and the pressure boundary). 

The pressure balance is determined for either the over or under expanded 
case with impact or ambient freestream pressure. 

CALLING SEQUENCE 

CALL HYPER (PB, I, K, ITYPE) 

where (PB) is the boundary pressure, (I, K) locates the boundary point, and 
(ITYPE) indicates if upper or lower boundary is involved. 

UTILITY ROUTINES AND COMMON REFERENCES 

COMMON/DATAR/ 

COMMON/STEPC/ 

FABLE 

POFEM 

EMOFV 

FNEWTN 

OVEREX 

ITSUB 

THETPM 

ERRORS 

METHOD OF SOLUTION 

The boundary pressure (may be impact or ambient) is compared to 
static pressure of the corner point. Depending on whether the comparison 
indicates the flow is over or under expanded, a branch is made to (OVEREX) 
or (THETPM). In either of these routines an iterative process balances the 
boundary pressure with the flow field pressure at the boundary. 
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SUBROUTINE NAME: ICHECK 
DESCRIPTION 

This function checks on the absolute difference in distance, Mach no. , 
and flow angle between the new point and both the right and left running base 
points. This check determines if a point is added, deleted, or left in as 
calculated. 

CALLING SEQUENCE: 

IDO = ICHECK (I, K, II, J, J2, K) 

where (I, K), (II, J), (12, K) designate the new characteristic point, the right 
running base point, and the left running base point, respectively. 

UTILITY ROUTINES AND COMMON REFERENCES 

common/dat ar/ 
common/divcri/ 

FABLE 

EMOFV 

VOFEM 

METHOD OF SOLUTION 


Not applicable. 
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SUBROUTINE NAME ; IDTAPE 
DESCRIPTION 

This subroutine writes the gas properties which were input via cards 
on the flow field program tape. The format used to write them on tape is 
compatible with that used for real gas. 

CALLING SEQUENCE 

CALL IDTAPE (UNITS) 

where (UNITS) indicates whether the gas properties are being read in with 
English (ENG) or MKS units. 

UTILITY ROUTINES AND COMMON REFERENCES 

COMMON/ DA tar/ 
common/ tapefo/ 
common/contrl/ 
common/head/ 

UTILITY - None 
METHOD OF SOLUTION 

Gas property data is read in from cards. If not already in MKS units, 
the data is converted. This converted data is then written on the. flow field 
tape (tape 13). 
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SUBROUTINE NAME: INITP 
DESCRIPTION 

This subroutine initializes the values of various control parameters, 
thereby providing for proper operation of the program. These initial values 
include: 

1. the counter for the upper and lower boundary equations, 

2. the counter for the first characteristic line, 

3. the initial number of degrees per Frandtl-Meyer ray, 

4. convergence criteria, 

5. maximum number of iterations. 

CALLING SEQUENCE 

CALL INITP 

UTILITY ROUTINES AND COMMON REFERENCES 

COMMON/CONTRL/ 

COMMON/CRITER/ 

common/datar/ 

common/head/ 

COMMON/ST EPC/ 

UTILITY - None 

METHOD OF SOLUTION 
Not applicable. 
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SUBROUTINE NAME: INRSCT 
DESCRIPTION 

This subroutine finds the intersection of two straight lines. 

CALLING SEQUENCE 

CALL INRSCT (Rl, XI, BETA1, R2, X2, BETA2, R3, X3) 

where (Rl, XI, BETA1) and (R2, X2, BETA2) define the equations of the two 
straignt lines which intersect at (R3, X3). 

UTILITY ROUTINES AND COMMON REFERENCES 
None 

METHOD OF SOLUTION 

The equations of the straight lines are written 

r = tan/?j(x - Xj) + r^ 

and 

x = cot $2 (r ~ r^) + 

These equations are solved for x but a test on the slopes is made to prevent 
indeterminate forms. If an indeterminate form is possible, the points are 
mapped one onto another, thus precluding the possibility of indeterminancy 
except when the lines are parallel. 
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FUNCTION NAME : ITERM 
DESCRIPTION 

ITERM tests each characteristic point to determine if it is within the 
predefined problem limits. If the point falls outside the limits, the present 
characteristic line is terminated. Control is then returned to PHASE1 for 
initiation of a new line. 

CALLING SEQUENCE 

FUNCTION = ITERM (IP, K) 

where (IP) identifies the characteristic point on the new (K) line. 

UTILITY ROUTINES AND COMMON REFERENCES 

COMMON/CUTFO/ 

COMMON/DATAR/ 

UTILITY - None 

METHOD OF SOLUTION 

The angular orientation of a line drawn from the upper or lower cutoff 
coordinates to the characteristic point is determined. Comparing this angle 
to the angle of the upper or lower cutoff line determines if the point is inside 
or outside the problem limits. 
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SUBROUTINE NAME; ITSUB 
DESCRIPTION 

This subroutine controls the iterative solution of any set of equations 
which can ultimately be expressed as a function of one variable. The routine 
can also be used to control an integration loop. 

CALLING SEQUENCE 

CALL ITSUB (FOFX, X. SAVE, CCNV, NTIMES) 

(FOFX) - function of X which is driven to zero 

(X) - variable which is iteratively solved for 

(SAVE) - program control, i.e., SAVE(1' control counter, 

SAVE(2) is X increment 

(CCNV) - convergence criteria 

(NTIMES) - maximum number of iterations 

UTILITY ROUTINES AND COMMON REFERENCES 
None 

METHOD OF SOLUTION 

ITSUB modifies (X) in the proper direction by the decrement value 
(SAVE(2) ) until the root has been bracketed. The method of false position 
is then used to modify X until the solution is reached. Immediately after 
entering ITSUB each time, the function is inspected for convergence. If 
the function has converged, a program control is set, and computer control 
is transferred to the calling routine. 
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SUBROUTINE NAME: JWBTSS 
DESCRIPTION 

This subroutine performs all the flow field tape search operations for 
determining flowfield properties at each start line point for the automatic re- 
start of reflected shock. The routine is entered with the coordinates of a 
point on the start line, and a linear interpolation with .n the characteristics 
mesh produces flow field data at the desired point. 

CALLING SEQUENCF 

CALL .1 • BTSS (POINT) 

where 

Point (1) = radial or vertical coordinate of desired point 

Point (2) = axial coordinate of desired point 

Point (3) = Mach number at desired point 

Point (4) = flow angle at desired point 

Point (5) = entropy level of flow at desired point 

Point (7) = velocity at desired point 

UTILITY ROUTINES AND COMMON REFERENCES 

COMMON/DUMPCO/ 

COMMON/FLAG/ 

READB 

READF 

TAPMOV 

METHOD OF SOLUTION 

The main function of the routine is to locate the start line point within 
the flow field or outside of its boundaries. If the start line point is found to 
be within the flow field, then flow field data is obtained at the location. This 
search operation is outlined in a series of steps performed by JWBTSS. 
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1. The routine first reads in the flow field boundary points and deter- 
mines if it is possible for the start line point to be inside the flow 
field. If this test is negative, a flag is set accordingly and control 
is returned to the calling routine J3LCKBX. 

If the test is positive, step 2 is taken. 

2. Characteristic lines are read in (one at a time) and retained in 
storage. As the lines are being read, each characteristic point is 
checked in an effort to bracket the body point. The search is con- 
ducted in a downstream direction, however, if the new point is up- 
stream of the previous point, the tape is backspaced six lines and 
bracketing is again attempted. If the tape is searched one hundred 
lines downstream of the starting point, or if the flow field boundary 
is reached, the tape is rewound and the search is started again. 

3. If the body point is bracketed, four to eight of the surrounding char- 
acteristic points are used in a linear interpolation to determine Mach 
number, flow angle, and entropy. 


Three flags are used as indicators of the results of the tape search. 


1. IF LAG 

2. JFLAG 

3. KFLAG 


1 

2 

1 

2 

1 

2 


indicates 

indicates 

indicates 

indicates 

indicates 

indicates 


boundary has not been read in. 

bound? ry has been read in. 

last start line point was not in flow field. 

last start line point was xn flow field. 

current start line point is not in flow field. 

current start line point is in flow field. 
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SUBROUTINE NAME: KIKOFF 
DESCRIPTION 

This subroutine allows control to return to the Main program if an 
error in the calculation is encountered, and controls the execution of mul- 
tiple cases. 

CALLING SEQUENCE 

CALL KIKOFF 

UTILITY ROUTINES AND COMMON REFERENCES 
COMMON/CONTRL/ 

common/quit/ 

DRIVER 

METHOD OF SOLUTION 
Not applicable. 
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SUBROUTINE NAME; LIMITS 
DESCRIPTION 

This subroutine tests the new boundary point to determine if it is within 
the limits of the current boundary equation. Depending on the test, the options 
are: 

1. use the current boundary equation, 

2. advance to the next boundary equation, or 

3. the current equation is the last one specified. 

CALLING SEQUENCE 

CALL LIMITS (I, K, ITYPE, IOK) 

where (I, K) represents the location of the boundary point in the PHO array, 
(ITYPE) tells if an upper or lower boundary is to be considered, and (IOK) 
is a control indicating if option 1, 2, or 3 is to be used. 

UTILITY ROUTINES AND COMMON REFERENCES 

common/ccntrl/ 

common/datar/ 

BOUND 

METHOD OF SOLUTION 

The radius (RMAX) and boundary angle (THETAMAX) at the limiting 
axial value (XMAX) is calculated in (BOUND). (RMAX) or (XMAX) is compared to 
(R) or (X) for the point in question. The results of the comparison determine 
if option 1, 2 or 3 is to be used. 


3-71 


LOCKHEED - HUNTSVILLE RESEARCH & ENGINEERING CENTER 


LMSC/HREC D 162220 -IV 

SUBROUTINE NAME: LIPIN 
DESCRIPTION 

LIPIN calculates information for the starting line points when the 
simplified straight start line option is used (i. e. , when ICON (2)^ 2). 

CALLING SEQUENCE 

CALL LIPIN (COOR, S, INTOT, DELM) 

where (COOR) is the starting line information array, (S) is the entropy 
level of the start line, (1_ TOT) is the total number of input points speci- 
fied (50 MAX), and (DELM) is Mach number gradient along the start line. 

UTILITY ROUTINES AND 3QMMON REFERENCES 

common/input/ 

COMMON/CONTRL/ 

RGVOFM 

UOFV 

METHOD OF SOLUTION 

The star* line input data is divided into the specified number of incre- 
ments. Radial gradients in Mach number, X, 0, are calculated. A circular 
arc transformation is applied to the input line data points to concentrate the 
points near the outer boundary. 
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SUBROUTINE NAME: MASCON 
DESCRIPTION 

This subroutine calculates the Mach number distribution at an area 
downstream of the throat such that total mass flow is conserved. Mass 
flow, calculated at the throat, is used as the constant for comparison. 

CALLING SEQUENCE 

CALL MASCON (E, SE, DELM) 

where (E) is the input line array (CORLIP), (SE) is the input line entropy 
level, and (DELM) is the Mach number gradient along the start line. 

UTILITY ROUTINES AND COMMON REFERENCES 

COMMON/CONTRL/ 

RGVOFM 

RHOFEM 

ERRORS 

EMOFV 

ITSUB 

METHOD OF SOLUTION 

The mass flow rate at the throat (fti*) is calculated. This (m*) is com- 
pared to that at the input line location for an initial Mach number distribution. 
The Mach number distribution is then perturbed until mass flow is conserved. 


3-73 


LOCKHEED - HUNTSVILLE RESEARCH & ENGINEERING CENTER 



LMSC/HREC D162220-IV 


SUBROUTINE NAME: MASSCK 
DESCRIPTION 

This subroutine keeps a running check on the mass flow. Mass flow 
at the starting line is calculated and compared with that crossing each 
characteristic line downstream. 

CALLING SEQUENCE 

CALL MASSCK (ILAST, ISTART, K) 

where (ILAST) is the last point on the characteristic line, (ISTART) is a 
counter for characteristic lines which emanate from the start line, and (K) 
represents the characteristic line under consideration. 

UTILITY ROUTINES AND COMMON REFERENCES 

COMMON/DAT AR / 

COMMON/INPUT / 

CCMMON/CON RL/ 

COMMON/ST EPC/ 

common/force/ 

common/wtsavf/ 

FABLE 

EMOFV 

RHOFEM 

UOFV 

METHOD OF SOLUTION 

The mass flow through the start line is calculated and stored. Mass 
flow through lines downstream is calculated and these values compared with 
the initial value. A percent change in mass flow is printed for each charac- 
teristic line. The total mass flow passing under each point on a characteristic 
line is stored so the mass flow can be written on the output tape to permit 
streamline tracing. 
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SUBROUTINE NAME ; MOCSOL 
DESCRIPTION 

This subroutine provides all two dimensional or axisymmetric method- 
of-characteristics solutions. The new point being solved for may be one of 
five possible types: 


1. 

interior point 

2. 

upper wall point 

3. 

upper free boundary point 

4. 

lower wall point 

5. 

lower free boundary point 

CALLING SEQUENCE 


CALL MOCSOL (I, K, II, Kl, 12, K2, IFLAG, ITYPE) 

where (I, K) identifies the storage location for the new point to be computed, 
(II, Kl]| identifies the right running (or upper boundary) known point, and 
(12, K2) identifies the left running (or lower boundary) known point. (IFLAG) 
is an error indicator and (ITYPE) selects the type calculation. 

UTILITY ROUTINES AND COMMON REFERENCES 


common/contrl/ 

BOUND 

common/criter/ 

ROTERM 

common/datar/ 

VOFEM 

common/gascon/ 

RGMOFP 

fable 

FNEWTN 

inrsct 

TOFV 

POFEM 

UOFV 

EMOFV 

RHOFEM 


ERRORS 


METHOD OF SOLUTION 

The four characteristic equations are written as a function of five 
variables (R,X,9,V,S). An additional relationship is obtained by assuming 
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the entropy (S) varies linearly between known data points. Using these 
characteristic equations in finite difference form, the routine solves for a 
new mesh point, knowing two mesh points of opposite family. 

The solution is begun by setting the average values of properties over 
the step length equal to the known values at the base points. Subsequent 
passes in the iterative solution result in •'updated" average values. The 
iterative solution is continued until the desired convergence on velocity or 
flow angle is reached or until the maximum number of iterations is exceeded. 

For a detailed derivation of the characteristic equations and a descrip- 
tion of their application in finite difference form to the solution of the char- 
acteristic mesh, see Vol.III, Section 6. 
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S UBROUTINE NAME; MONO 
DESCRIPTION 

This subroutine determines if the characteristic point locations are 
monotonic along left or right running lines. 

CALLING SEQUENCE 

CALL MONO (II, Jl, 12, J2, 13, K3, IOK) 

where (II, Jl), (12, J2), (13, K3) designate the right running base point, the 
left running base point, and the new characteristic point, respectively. ICK 
is a flag returned to the calling program which is (0) if solution is monotonic, 
(1) if a non-monotonic condition has occurred along the right running charac- 
teristic, and (2) if one has occurred along a left running line. 

UTILITY ROUTINES AND COMMON REFERENCES 

common/da tar/ 

UTILITY - None 
METHOD OF SOLUTION 

Envelope shock waves are detected in a method-of-characteristics 
solution by crossing of characteristic lines which, mathematically, causes 
a discontinuity in the flow properties. This discontinuity is interpreted 
as a shock wave. This routine is supplied the two base points and the 
resultant new mesh point. Using this information, a discontinuity, if it 
exists, is detected. 
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Now 

^ fTT ~ ^ix ^1 ” ^x) ^ ~ Mj) j| 

and 

*m = -n + *u { c ° s < 5 n + 5 n )‘ + sin ‘®n < m u ) j) 

Effectively, this routine checks the sign of X^, ^ and sets IOK accordingly. 


if 

*r a h 

c. 0 

IOK = 

0 

if 

>* 

hi 

< 0 

IOK = 

1 

if 


< 0 

IOK = 

2 
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SUBROUTINE NAME : OUT 
DESCRIPTION 

This subroutine writes the calculated data at characteristic points 
along with the corresponding title and headings. 

CALLING SEQUENCE 

CALL OUT (II, 12, K) 

where (II, 12) ref *r to the point numbers of the points to be output (any 
number of points may be output at one time)* (K) represents the current 
characteristic line (takes on the value 1 or 2)* 


UTLuITY ROUTINES AND COMMON REFERENCES 


common/con TRL,/ 

EMOFV 

common/dat ar/ 

TJOFF.M 

COMMON/GASCON/ 

POFEM 

common/head/ 

RHOFLM 

PAGE 

TOFEM 

FABLE 

VOFEM 


ESHOCK 


METHOD OF SOLUM ION 
Not applicable. 
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SUBROUTINE NAME: OUTBIN 
DESCRIPTION 

This subroutine writes the calculated characteristic point data on th« 
binary output tape. This is done for any number of characteristic points. 

CALLING SEQUENCE 

CALL OUTBIN (II, 12, J) 

where (II, 12) identifies the range of points to be written on tape (II is fir 
point, 12 is last). (J) represents the current characteristic line (1 or 2). 

UTILITY ROUTINES AND COMMON REFERENCES 

COMMON/D AT AR / 

COMMON/WTSAVF/ 

COMMON/ITAPE/ 

FABLE 

EMOFV 

METHOD OF SOLUTION 
Not applicable. 
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SUBROUTINE NAME: OVEREX 
DESCRIPTION 

This subroutine solves for the shock angle at the nozzle lip when the 
flow is over expanded. Provisions are made to calculate the shock angle 
for an upper or lower lip point. Real gas effects are considered in calcu- 
lating flow properties downstream of the shock. 

CALLING SEQUENCE 

CALL OVEREX (PB, I, K, ITYPE) 

where (PB) is the freestream pressure at the boundary, (I, K) defines the 
location of the lip point in the characteristic data (PHO) array and (ITYPE) 
indicates whether an upper or lower boundary is to be considered. 

UTILITY ROUTINES AND COMMON REFERENCES 

COMMON/DAT AR/ 

FABLE 

EMOFV 

ESHOCK 

POFEM 

ITSUB 

UOFV 

METHOD OF SOLUTION 

An initial shock angle is assumed. This shock angle is perturbed in 
ITSUB and the result used to calculate flow properties downstream of the 
shock, including static pressure. The calculated static pressure is com- 
pared with the boundary pressure and the difference is driven sufficiently 
close to zero to satisfy convergence criteria. 
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SUBROUTINE NAME: PAGE 
DESCRIPTION 

This subroutine page ejects and writes the header comments and page 
number on each page of printout. 

CALLING SEQUENCE 

CALL PAGE (LCNT) 

where (LCNT) is a counter which monitors the number of lines of printed 
output per page. (LCNT) is re -initialized in PAGE. 

UTILITY ROUTINES AND COMMON REFERENCES 

common/head/ 

common/contrl/ 

UTILITY - None 

METHOD OF SOLUTION 

When the maximum number of lines per page (55) have been* output, (PAGE) 
is called to page eject. It then prints the identifying information and the page 
number, increments the page number and re-initializes the line counter. 
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SUBROUTINE NAME: PHASE1 


DESCRIPTION 

This subroutine provides the necessary logic to successively employ 
the proper calculation at the proper time in order to describe the entire 
characteristic mesh. Direction is given to calculate the flow field through- 
out the nozzle and plume, and termination is achieved when a right running 
shock intersects the boundary or problem limits have been reached. 

CALLING SEQUENCE 

CALL PHASE 1 (IF INIS) 

where (IFINIS) is a flag to bring in additional logic routines for restarting 
the characteristic solution for reflecting a shock which has intersected the 
axis of symmetry. 


UTILITY ROUTINES AND COMMON REFERENCES 


COMMON/CONTRL/ 

UOFV 

COMMON/GASCON/ 

ITERM 

COMMON/DATAR/ 

RGVOFM 

COMMON/ITAPE/ 

HYPER 

C ommon/in pu t/ 

POFEM 

common/del/ 

EMOFV 

common/stepc/ 

RGMOFP 

common/which/ 

VOFEM 

common/finite/ 

THETPM 

FABLE 

MONO 

OUT 

TURN 

THRUST 

OVEREX 

OUTBIN 

ERRORS 

MOCSOL 

SHOCK 

LIMITS 

MASSCK 

BOUND 

ICHECK 

PR AND T 



ESHOCK 

SAVER 
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METHOD OF SOLUTION 

In general terms, the flow properties at an unknown characteristic 
point may be expressed as some operation (ij>) on some number • • • » 

*m> of known points* These operations (0) differ according to the type of 
unknown point to be calculated. Presently, the six types of unknown points 
dealt with or six operations ( 0) performed are: 

1. starting line point (0q) 

2. interior point (0^) 

3. boundary point (upper or lower) (0g) 

4. attached shock point (0^g) 

5. shock point (0g) 

6. Prandtl- Meyer point (upper or lower, solid or free) (0pi v i) 

Basically, PHASE 1 contains the fixed point arithmetic necessary to perform 
the above mentioned calculations. Through this system of fixed point arith- 
metic, the necessary types of calculations are performed on the proper known 
characteristic points to produce the new characteristic point. For a detailed 
description of the fixed point arithmetic and examples of its use, see 
Appendix. 
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SUBROUTINE NAME: PLUMIN 
DESCRIPTION 

This subroutine reads in the input data (input via cards) necessary to 
perform the method-of-characteristics solution. This input information is 
routed by PLUMIN to various supporting routines depending on the options 
selected. 

CALLING SEQUENCE 

CALL PLUMIN 

UTILITY ROUTINES AND COMMON REFERENCES 


COMMON/CONT RL/ 

GASRD 

VOFEM 

common/cutfo/ 

BOUND 

RGMOFB 

common/datar/ 

LIPIN 

CHECK 

common/gas con/ 

AOASTR 

ITSUB 

common/head/ 

MASCON 

ERRORS 

common/input/ 

RGVOFM 

VISCUS 

COMMON/ST EPC 

UOFV 


COMMON/POINTS/ 

POFEM 


COMMON/DIVCRI/ 

PLMOUT 


common/quit/ 

THROAT 


common/which/ 

SOLVE 



METHOD OF SOLUTION 
Not applicable. 


3-85 


LOCKHEED - HUNTSVILLE RESEARCH & ENGINEERING CENTER 



LMSC/HREC D162220-IV 


SUBROUTINE NAME: PLMOUT 
DESCRIPTION 

This subroutine prints the data read by (PLUMIN). 

CALLING SEQUENCE 

CALL PLMOUT (KP, LCNT) 

where (KP) is a control parameter which is set in PLUMIN, and (LCNT) is 
the printed line counter. 

UTILITY ROUTINES AND COMMON REFERENCES 

COMMON/CONTRL/ 

COMMON/CUTFO/ 

COMMON/TAPEFO/ 

common/datar/ 

COMMON/GASCON/ 

COMMON/HEAD/ 

common/input/ 

common/divcri/ 

PAGE 

FABLE 

EMOFV 

METHOD OF SOLUTION 
Not applicable. 
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FUNCTION NAME ; POFEM 
DESCRIPTION 

This function computes the local static pressure as a function of Mach 
number and entropy. 

CALLING SEQUENCE 

P = POFEM (EM, S) 

where (P) is the resultant static pressure found from the Mach number (EM) 
and entropy (S). NOTE : The appropriate values of the gas properties must 
be stored in common upon entry to this routine. 

UTILITY ROUTINES AND COMMON REFERENCES 

common/gascon/ 

UTILITY - None 
METHOD OF SOLUTION 

Thermally perfect gas relationships are used to find the pressure. 

- y 

^-s/r /, . 7 - 1 „,2\ y_1 
p = p 0 6 l 1 + —i — M ) 
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SUBROUTINE NAME : PRANDT 
DESCRIPTION 

This subroutine computes the Prandtl-Meyer expansion angle for a 
given boundary angle and divides this angle into a series of expansion 
"rays” (unless the number of rays has been specified in the input). The 
flow properties at each angular increment are set and stored in the Char- 
acteristic Data (PHO) array. 

CALLING SEQUENCE 

CALL PRANDT (I, J, THETAB, NPM, IF LAG, ITYPE) 

where 

(I) - represents the corner point 

(J) - indicates a characteristic line 
(THETAB)- is the boundary angle 

(NPM)- number of Prandtl-Meyer increments (calculated in PRANDT) 
(IF LAG)- error flag 

(ITYPE)- indicates if upper (2) or lower (1) boundary 

UTILITY ROUTINES AND COMMON REFERENCES 

COMMON/CRIT ER/ 

COMMON/DAT AR/ 

COMMON/GASCON/ 

COMMON/STEPC/ 

COMMON/CONTRL/ 

FABLE 

THETPM 

UOFV 

METHOD OF SOLUTION 

The routine is entered with known flow properties at the point of dis- 
continuity along with the known corner and boundary flow angles. From the 
known angles and the preset number of degrees per ray, the number of 
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increments is calculated. The distribution of P-M rays is then adjusted by 
a weighting function. Subroutine (THETPM) is entered with known initial 
conditions and the number of degrees per ray and returns with a velocity. 
These new conditions are then set into the (PHO) array. 
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SUBROUTINE NAME: READB 
DESCRIPTION 

Subroutine READB reads the position of the boundary points of the flow 

field. 


This information is read from the binary output tape for setting up the 
start line for automatically restarting a reflected shock. 

Subroutine JWBTSS uses this to determine boundary position of flow 

field. 

CALLING SEQUENCE 

CALL READB (X,R,ITOTl) 

X = X coordinate of last point on characteristic line 
R = R coordinate of last point on characteristic line 
ITOT1 = Number of points on the line 

UTILIY ROUTINES AND COMMON REFERENCES 
None 

METHOD OF SOLUTION 
Not applicable 
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SUBROUTINE NAME : READF 
DESCRIPTION 

SUBROUTINE READF reads one characteristic line from flow field tape 
and saves the following data at each point: 

1. X position of the point, 

2. R position of the point, 

3. Mach number at the point, 

4. Flow angle at the point, 

5. Entropy at the point, 

6. Velocity at the point. 

This information is read from a binary output tape for setting up the 
start line for the automatic restart of a reflected shock. 

CALLING SEQUENCE 

CALL READF (J, ITOT1) 

J — number of characteristic line 

ITOT1 — number of points in characteristic line 

UTILITY ROUTINES AND COMMON REFERENCES 

common/dumpco/ 

UTILITY - None 

METHOD OF SOLUTION 
Not applicable 
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FUNCTION NAME : RGMOFP 
DESCRIPTION 

This subroutine finds Mach number as a function of pressure, O/F ratio 
and entropy. The difference between this routine and EMOFP is that in this 
case the gas properties are not known prior to entry. 

CALLING SEQUENCE 

EM = RG\iOFP (P, OF, S) 

where (EM) is the resultant Mach number, (P) is the local static pressure, 
while (S) is the local entropy^ and (OF) is the local O/F ratio. 


UTILITY ROUTINES AND COMMON REFERENCES 


common/contrl/ 

common/gascon/ 

POFEM 

EMOFV 

ITSUB 

FABLE 


COMMON./ DATAR/ 

common/tapefo/ 

VOFEM 

EMOFP 

ERRORS 


METHOD OF SOLUTION 

The real gas tables have, as independent variables, OF ratio, entropy 
and velocity. If the velocity is not known, an iterative solution must be 
employed to find Mach number from pressure, entropy, and OF ratio. 
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FUNCTION NAME : RGVOFM 
DESCRIPTION 

This subroutine finds velocity as a function of Mach number, entropy 
and O/F ratio. The difference between this routine and VOFEM is that the 
gas properties are not known prior to entry. 

CALLING SEQUENCE 

Y = RGVOFM (OF, S, EM) 

where V is the resultant velocity computed from OF ratio (OF), entropy (S), 
and Mach number (EM). 

UTILITY ROUTINES AND COMMON REFERENCES 


FABLE 

common/contrl/ 

VOFEM 

common/tapefo/ 

EMOFV 

common/datar/ 

ITSUB 


ERRORS 


METHOD OF SOLUTION 



The real gas tables have, as independent variables, OF ratio, entropy 
and velocity. If the velocity is not known, an iterative solution must be 
employed to find the velocity from Mach number, OF ratio and entropy. 
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FUNCTION NAME: RHOFEM 


DESCRIPTION 

This function computes the local density as a function of Mach number 
and entropy. 

CALLING SEQUENCE 

RHO * RHOFEM (EM, S) 

where RHO is the resultant density found from local Mach number and local 
entropy. NOTE: The appropriate values of the gas properties must 

be stored in common upon entry to this routine. 

UTILITY ROUTINES AND COMMON REFERENCES 

common/gascon/ 

POFEM 

METHOD OF SOLUTION 

Thermally perfect gas relationships are used to find the density. 

p * Mi + 
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FUNCTION NAME : ROTERM 
DESCRIPTION 

This routine computes the geometrical factor (Fj, F II> used in the axisym- 
metric term of the capability equation and as an interpolation parameter. 

CALLING SEQUENCE 

F = ROTERM (THETA, DELTA, EMU, R3, RI) 

where (THETA) is or 9 jj (flow angles of the known points) 

(DELTA) selects quadrant 

(EMU) is JIj or JJjj (Mach angles of the known points) 

(R3) is "Fjjj or Xjjj (coordinates of new point) 

(RI) is rj or X j (coordinates of known point) 

UTILITY ROUTINES AND COMMON REFERENCES 
None 

METHOD OF SOLUTION 

The method-of-characteristics solution uses this routine to determine 
a coefficient needed in its solution. This term (see Equation (6-29), Section 
6 of Reference 1) can be written as: 

F = I 8 H < d IM - d > 

+ «(9 + Ji - f)) 

By the proper choice of d (r or x), 6 and the sign of , indeterminant 
forms are eliminated in the evaluation. 
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SUBROUTINE NAME: SAVER 
DESCRIPTION 

This routine saves the flow properties at the nozzle lip for setting up 
the start line for the automatic restart of the reflected shock case. 

CALLING SEQUL-.CE 

CALL SAVER (I, K) 

where I is the lip point number on the characteristic line, and 

K is the characteristic line on which the lip point was found. 

UTILITY ROUTINES AND COMMON REFERENCES 

common/dat ar/ 
common/shckpt/ 

UTILITY - None 

METHOD OF SOLUTION 
Not applicable 
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SUBROUTINE NAME; SETUP 
DESCRIPTION 

For the automatic restart of a reflected shock this subroutine locates 
and distributes the start line points along a straight line from the nozzle lip 
or first point on the original start line to the point where the shock has inter- 
cepted the axis. 

CALLING SEQUENCE 

CALL SETUP (COOR, INTOT) 

where COOR is the starting line information array, and 
INTOT is the number of start line points specified. 

UTILITY ROUTINES AND COMMON REFERENCES 

common/input/ 

common/'contrl/ 

UTILITY - None 
METHOD OF SOLUTION 

The reflected shock start line data is divided into the specified number 
of increments. Radial gradients in R, X and 0 are calculted. A circular arc 
transformation is applied to the input line data to concentrate the points near 
the outer boundary. 
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FUNCTION NAME: SHOCK 


DESCRIPTION 

This subroutine iteratively adjusts the shock strength in order to 
satisfy the oblique shock relations and the flow field properties simultaneously. 
The six different options are: 


1. interior right running shock wave 

2. interior left running shock wave 

3. right running shock wave at wall 

4. left running shock wave at wall 

5. right running shock reflected from boundary 

6. left running shock reflected from boundary 
NOTE: Options 2, 4 and 6 currently not executed. 

CALLING SEQUENCE 

CALL SHOCK (IN, KN, IKL, INI, JN, UH, I8FIN, I7FIN, IF LAG, ITYPE.IOPT) 


where 


(IN, KN) 
IKL 

(INI, JN) 
IJH 
18 FIN 
17 FIN 
IF LAG 
ITYPE 


is the location of the virtual point 

is the first point on the KN array 

is the location of the known shock point 

is the last point on the JN array 

is the final location of the upstream shock point 

is the final location of the base point on the downstream side 

is an error flag 

selects the type of calculation 


ITYPE 


IOPT 


1 1 for case (1) 

12 for case (2) 

21 for case (3) 

22 for case (4) 

31 for case (5) 

32 for case (6) 

Parameter for controlling direction and magnitude of 
perturbations during iteration. 
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UTILITY ROUTINES AND COMMON REFERENCES 


common/ contrl/ 
common/crit er/ 
common/gascon/ 
common/da tar/ 
common/irfec/ 

INRSCT 

FABLE 


ERRORS 

UOFV 

MOSCOL 

ROTERM 

ITSUB 

ESHOCK 


METHOD OF SOLUTION 

This subroutine is used to calculate the properties across an oblique 
shock wave as a function of the local characteristic lattice. In its operation, 
this routine instructs the calling routine as to the necessity of adjusting the 
counting scheme for network construction due to the presence of the shock 
wave. Diagrams for the six options are given below. 



( 2 ) 



Case (2), ITYPE = (12) 
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SUBROUTINE NAME; SOLVE 
DESCRIPTION 

This subroutine solves for the coefficients of the nozzle wall cubic 
equation described by points input to the program. 

CALLING SEQUENCE 

CALL SOLVE (ITYPE, M. YP) 

where ITYPE determines if the solution is for an upper or lower wall, 

M is the total number of input points describing the wall, and 

YP is the slope at the first point describing the wall. 

UTILITY ROUTINES AND COMMON REFERENCES 

A series of third order curve fits is performed. Each curve fit uses 
3 consecutive input points and the slope at the first input point. The coefficients 
for a polynomial wall equation are then determined from each curve fit. 

METHOD OF SOLUTION 


Not applicable. 
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SUBROUTINE NAME; SWITCH 
DESCRIPTION 

This routine transfers the upper boundary equations into the lower 
boundary array and the lower boundary equations into the upper boundary 
array. The appropriate signs on the equation coefficients are also changed 
so that the problem can be automatically restarted for a reflected shock case. 

calling sequence 

CALL SWITCH 

UTILITY ROUTINES AND COMMON REFERENCES 

COMMON/CONTRL/ 
common/dat ar/ 

UTILITY - None 

METHOD OF SOLUTION 
Not applicable. 
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SUBROUTINE NAME: TABLE 
DESCRIPTION 

This subroutine utilizes real or ideal gas information obtained from a 
master tape or input cards to calculate properties locally in the flow. The 
maximum size of the array used by (TABLE) is limited to five gas properties 
(V, R, y, T q , F ) at thirteen velocity "cuts" for each of two entropy cuts 
and 10 O/F cuts. 

CALLING SEQUENCE 

CALL TABLE (SS, VV, IF) 

where (SS) is the local entropy, (IF) is the O/F table of interest and (VV) is 
the local velocity at the point of interest. 

UTILITY ROUTINES AND COMMON REFERENCES 

common/xsicom/ 
common/contrl/ 
common/dat ar/ 
common/gascon/ 
common/fab/ 

TOFV 

POFEM 

EMOFV 

METHOD OF SOLUTION 

The routine is entered with an o/F table (IF), the local entropy (SS) and 
velocity (VV). A test is then made to determine if the gas is real or ideal. If 
the test indicates an ideal gas, the local properties are set to those stored in 
the (TAB) common array. If the test indicates real gas, a double interpolation 
scheme is utilized to locate gas properties between tabulated values of velocity 
and entropy. In the case of an entry beyond the range of the tables, an ideal 
gas extrapolation from the last table value is made to determine the gas 
properties. 
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SUBROUTINE NAME : TAPMOV 
DESCRIPTION : 

This subroutine moves the binary flow field tape through the gas property 
data to the start of the flow field information. 

CALLING SEQUENCE 

CALL TAPMOV 

UTILITY ROUTINES AND COMMON REFERENCES 

common/irecd/ 

UTILITY - None 

METHOD OF SOLUTION 
Not applicable. 
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SUBROUTINE NAME ; THETPM 
DESCRIPTION 

This subroutine performs a numerical integration to calculate 
properties through a Prandtl-Meyer expansion. Either the case of known 
final velocity or known final expansion angle may be handled. 

CALLING SEQUENCE 

CALL THETPM (OF, S, DELTA, VF, VI, IT, ITYPE) 

where (OF) is the local OF ratio 

(S) is the local entropy level 

(DELTA) is the total expansion angle 

(VF) is the final velocity downstream of the expansion 

(VI) is the initial velocity upstream of the expansion 

(IT) is a control parameter indicating if expansion to a solid 

wall or free boundary is taking place 

(ITYPE) indicates if upper (2) or lower (1) boundary 

UTILITY ROUTINES AND COMMON REFERENCES 

common/gascon/ ITSUB 

common/s tepc/ TOFV 

FABLE ERRORS 

METHOD OF SOLUTION 

The integral equation 

V F 

f ^M 2 - i ^ - ae = f(Vp) = o 
v i 

(M 2 = V 2 /VRT) 
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is solved knowing either the final velocity (V^) or the expansion angle (A0). 

As can be seen, if the final velocity (Vjr*) is known, the integration progresses 
straightforwardly to a solution. However, if the expansion angle is known, 
an iterative procedure must be employed to pick the velocity which produces 
the desired expansion angle. 
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SUBROUTINE NAME ; THROAT 
DESCRIPTION 

This subroutine computes the nozzle throat equation coefficients. 

CALLING SEQUENCE 

CALL THROAT (ITYPE, YP) 

where (ITYPE) is a 1 for an upper wall throat and a 2 for a lower wall 
throat. (YP) is the slope of the first point on the next curve. 

UTILITY ROUTINES AND COMMON REFERENCES 
COMMON/DAT ar/ 

common/points/ 

UTILITY - None 
METHOD OF SOLUTION 

The nozzle throat radius, radius of curvature, divergence angle and 
an axial corrdinate shift are used to find the coefficients for the conic 
equation: 

R = A {Vb + CX + 2DX 2 + E 
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SUBROUTINE NAME; THRUST 
DESCRIPTION 

This routine computes the vacuum thrust produced by a two-dimensional 
or axisymmetric nozzle. Addition of the thrust at the throat and the integrated 
pressure along the nozzle wall yields the final thrust, 

CALLING SEQUENCE 

CALL THRUST (I, K, II, Jl, ITYPE, ICALC) 

where (I, K) designates the unknown characteristic point and (II, Jl) is the 
known characteristic point. (ITYPE) specifies if the point is on the upper 
or lower boundary and (ICALC) is a counter with the values of 1, 2 or 3. 

(1 specifies integration at the throat, 2 - along the nozzle and 3 - at the exit.) 

UTILITY ROUTINES AND COMMON REFERENCES 

common/contrl/ 

common/da tar/ 

common/force/ 

common/input/ 

common/stepc/ 

fable 

emofv 

POFEM 

RKOFEM 

METHOD OF SOLUTION 

Thrust is found by first computing the momentum thrust in the sonic 
area or throat of the nozzle. The static pressure is then integrated along the 
nozzle wall and the total thrust found by summation of the pressure and 
momentum terms. Inclusion of a factor in the incremental force term 
accounts for either two-dimensional or axisymmetric flow. 
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FUNCTION NAME : TOP EM 
DESCRIPTION 

This function computes the local static temperature as a function of 
Mach number. The gas properties at the point of interest are known prior 
to entry. TOFEM and TOFV are quite similar; the difference being if Mach 
number or velocity is the known quantity. 

CALLING SEQUENCE 

T = TOFEM (EM) 

where (T) is the one-dimens ionally calculated local static pressure which 
exists at the Mach number (EM). 

UTILITY ROUTINES AND COMMON REFERENCES 

common/gascoiV 

METHOD OF SOLUTION 

The calorically perfect gas relationship 



is solved for static temperature at the local Mach number. 
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FUNCTION NAME : TOFV 
DESCRIPTION 

This function computes the local static temperature as a function of 
velocity. The gas properties at the point of interest are known prior to entry, 
TOFV and TOFEM are quite similar; the difference being if Mach number or 
velocity is the known variable. 

CALLING SEQUENCE 

T = TOFV (V) 

where (T) is the one -dimensionally calculated local static pressure which 
exists at the velocity (V). 

UTILITY ROUTINES AND COMMON REFERENCES 

common/gascon/ 

UTILITY - None 

METHOD OF SOLUTION 

The calorically perfect gas relationship 



is solved for static temperature at the local velocity. 
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SUBROUTINE NAME: TURN 
DESCRIPTION 

This subroutine solves for a shock wave which has a known turning 
angle ( <J ). A condition of known turning angle exists when the flow is turned 
through a compression corner on a solid boundary. Real gas effects are 
considered in calculating conditions downstream of the shock. 

CALLING SEQUENCE 

CALL TURN (PU, PD, DELTA, IFLAG) 

where (PU, PD) represent flow conditions upstream and downstream of the 
shock, (DELTA) is the turning angle, and (IFLAG) indicates if the solution 
will converge or not. 

UTILITY ROUTINES AND COMMON REFERENCES 

common/criter/ ESHOCK 

common/contrl/ itsub 

COMMON/FINITE/ CONE 

FABLE UOFV 

EMOFV 
UOFEM 

METHOD OF SOLUTION 

An initial shock angle is guessed. This shock angle is used to calcu- 
late a turning angle. The calculated turning angle is compared to the known 
turning angle and successive iterations on shock angle are performed until 
the turning angle difference is sufficiently close to zero. 
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FU NCTION NAME : UOFEM 
DESCRIPTION 

This function computes the Mach angle at a local Mach number. A 
test is made to ensure that the Mach number is greater than one; prior to the 
calculation. 

CALLING SEQUENCE 

EMU = UOFEM (EM) 

where (EMU) is the Mach angle which exists at the local Mach number (EM). 

UTILITY ROUTINES AND COMMON REFERENCES 

ERRORS COMMON - None 

KIKOFF 

METHOD OF SOLUTION 

The following equation 



is solved for the local Mach angle. 
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FUNCTION NAME: UOFV 

l , , i — 1 1,, . i. 


D ESCRIPTION 


This function computes the Mach angle at a local velocity. 


• CALLING SEQUENCE 

EMU = UOFV (V) 

where (EMU) is the Mach angle which exists at the local velocity (V). 

UTILITY ROUTINES AND COMMON REFERENCES 
COMMON - None 
UOFEM 
EMOFV 


METHOD OF SOLUTION 


The local velocity is converted into a Mach number using (EMOFV). 
Function (UOFEM) is then entered with the calculated Mach number and the 
Mach angle obtained from the following equation. 


p = tan 




1 


V 


M 2 - 1 
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SUBROUTINE NAME; VISCUS 
DESCRIPTION 

This subroutine calculates the laminar or turbulent boundary layer 
thickness and velocity distribution at the nozzle exit. This boundary layer 
adjustment is made on the exit plane starting line. Only the supersonic 
portion of the boundary layer is considered. 

CALLING SEQUENCE 

CALL VISCUS (N. NBL, XL, CU) 

where (N) is the power of the velocity profile (v/v e( jg e = (y/y e< jg e )*^) . 
(NBL) is number of boundary layer points specified. (XL) is a characteristic 
length (usually the nozzle length). (CU) is a conversion factor for mixed 
units of length in the boundary layer calculation. 

UTILITY ROUTINES AND COMMON REFERENCES 


common/input/ 

EMOFV 

common/ contrl/ 

TOFEM 

common/gascon/ 

POFEM 

fable 

VOFEM 


UOFV 

METHOD OF SOLUTION 



The number of starting line points within the supersonic portion of the 
boundary layer is specified. The velocity profile for a laminar or turbulent 
boundary layer is calculated and distributed to the starting line points. These 
points are then transferred to a right running characteristic line and the 
remaining inviscid portion of the starting line is attached. 
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FUNCTION NAME ; VOFEM 
DESCRIPTION 

This function computes velocity as a function of Mach numhjer. Ideal 
gas relations are used and the gas properties are known prior to entry. 

CALLING SEQUENCE 

V s VOFEM (EM) 

where (V) is the local velocity which corresponds to the local Mach number 
(EM). 

UTILITY ROUTINES AND COMMON REFERENCES 

COMMON/GASCON/ 

TOFEM 

METHOD OF SOLUTION 

The ideal gas relationship 

/Ry (T - T) 

is solved for velocity. Local static temperature (T) is obtained from the 
input Mach number. 
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SUBROUTINE NAME; WEAK 
DESCRIPTION 

This subroutine determines the independent variables (SD, VD) down- 
stream of a weak oblique shock. The gas properties upstream of the shock 
are known prior to entry. 

CALLING SEQUENCE 

CALL WEAK (OF, SU, VU, ERS, DELTA, SD, VD) 

where (OF) is the upstream O/f ratio, (SU, VU) are the upstream entropy 
and velocity, (EPS, DELTA) are the shock angle and turning angle, and (SD, 

VD) are the downstream entropy and velocity. 

UTILITY ROUTINES AND COMMON REFERENCES 

COMMON/GASCON/ 

TABLE 

EMOFV 

POFEM 

RHOFEM 

ENTROP 

DELTAF 

METHOD OF SOLUTION 

From the known upstream entropy and velocity, the local gas properties, 
pressure, density, and upstream Mach number are calculated. The entropy 
rise across the shock is added to the upstream entropy to get total downstream 
entropy. Downstream velocity is calculated from the following relationship. 

V cos (0 

V = 

D cos (£ - 6) 
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FUNCTION NAME : WOFA 
DESCRIPTION 

This function computes the weight flow per unit area as a function of 
Mach number. This calculation is only used in function AOASTR. 

calling sequence 

WEIGHT FLOW « WOFA (EM) 
where (EM) is the local Mach number. 

UTILITY ROUTINES AND COMMON REFERENCES 

common/gascon/ 

UTILITY - None 
METHOD OF SOLUTION 

Weight flow per unit area (w/A) is calculated from ideal gas relations. 
The equation used is 
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3.5 PROGRAM OUTPUT 


This subsection contains a description of the output scheme utilized 
by the program and a description of the error messages printed out by the 
program. 


3-118 


LOCKHEED - HUNTSVILLE RESEARCH & ENGINEERING CENTER 



p 


LMSC/HREC DI62220-IV 


3.5.1 Description of Program Output 

The methods of characteristics program output is organized in a logi- 
cal fashion with the data presented in an easily understood form. The initial 
pages consist of a printout of the input data including the real gas tables ob- 
tained from the master tape. Characteristic data are organized along left- 
running characteristic lines with all the pertinent information printed for 
each characteristic point on each characteristic line. Numbered flags on 
the example printout sheets (pp. 3-124 through 3-127) correspond to numbered 
comments listed below. 

GROUP I - IDENTIFICATION 

(T) Case Number: Appears on each page - may be a maximum of five digits. 

( 2 ) Title: Identifies particular run, appears on each page and may be 72 
spaces. 

GROUP 2 - RUN CONTROL 

(J) Run Control Parameters: These 16 parameters control the execution 
of the program according to the options selected. (See Input/Output 
FORTRAN symbols section for explanation of individual parameters.) 

GROUP 3 - BOUNDARY EQUATIONS (See Input/Output FORTRAN Symbols 
for detailed description.) 

^4) Tvpe Equation: Identifies type of boundary equation selected. (1 - conic, 
2 - poly, 3 - free bound) 

® ITRANS : Indicates if a discontinuity follows this equation, 

0 - no discontinuity, 1 - discontinuity). 

(6) Coefficients : Apply to upper or lower boundary equations. 

(V) XMAX : Maximum value of (x) for which present equation applies. 
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GROUP 4 - GAS IDENTIFICATION 


8) Gas Identification : Identifies gas on master tape for which gas table 
is printed. 


GROUP 5 - REAL/lDEAL GAS PROPERTIES (See Input/Output FORTRAN 
Symbol section) 


0 

© 

© 

(l2) 


© 


O/f Ratio : Th:* O/F ratio for the particular table. 

Entropy Cuts ; May be two maximum for each O/F table, value is 
relative to first entropy level in each table. 

Velocity Cuts : May be 13 maximum at each entropy cut. 

M Gas Constant 11 ; Value associated with the particular velocity, entropy 
and O/F ratio. 

Isentr^pic Exponent : Value associated with the particular velocity, 
entropy and O/F ratio. 

Static Temperature : Value associated with the particular velocity, 
entropy and O/F ratio. 

Static Pressure : Value associated with the particular velocity, entropy 
and O/F ratio. 

Note : Ideal gas format is similar (O/F) is as read in from card, (S) 
and (V) are generally zero and only one value of R, y, T , P q 
is printed. 


GROUP 6 - STARTING LINE INFORMATION (See Input/Output FORTRAN 
Symbol Section) 

@ R ; Radial distance to characteristic point (may be any unit of length 

or non-dimensional). 

(T7) X: Axial distance to characteristic point (may be any unit of length 

or non-dimensional). 

Note : The only restriction on units for X and R is that they be consistent, 
however thrust and mass flow calculations assume feet. 
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@ Mach Number : Local value at particular R and X location (may be any 
supersonic value). 

(T 9 ) Flov/ Angle : Local value at particular R and X location (degrees). 

^ 2 2 
(20) Entrcpy : Local value at particular R and X location (ft /sec -°R). 

(TT) O/f Ratio ; Local value at particular R and X location. 

Note : The starting line may be obtained using various options (see 
Input Guide); however, the format on the printout remains the 
same. 

GROUP 7 - MESH CONTROL CRITERIA (See Input/Output FORTRAN Symbol 
Section DIV(l) - DIV(6)) 

(22) DLI : Length insertion criteria. 

DLL : Length deletion criteria. 

DMI : Mach number insertion criteria. 

DML : Mach number deletion criteria. 

DTI : Flow angle insertion criteria. 

DTL : Flow angle deletion criteria. 

' 8 - RUN CUTOFF INFORMATION (See Input/Output FORTRAN Symbol 
Section CUTDAT(l) — CUTDAT(6) 

II: Radial coordinate of upper cutoff (units same as R on starting line). 
X; Axial coordinate of upper cutoff (units same as X on starting line). 
THETA: Angle of upper cutoff line (degrees) 

R_: Radial coordinate of lower cutoff. 

X: Axial coordinate ol’ lower cutoff. 

THETA: Angle of lower cutoff line. 
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GROUP 9 - TYPICAL LEFT CHARACTERISTIC LINE 





Characteristic Line: Identifies line for which data is printed. 
Characteristic Point : Identifies point for which data is printed. 

De scription : Describes the type of point. The options are: 

a. input point 

b. interior point (prints blank) 

c. wall point 

d. free boundary point 

e. shock point 

f. Prandtl -Meyer point. 

]R: Radial distance to the characteristic point. 

X: Axial distance to the characteristic point. 

M: Mach number at the characteristic point. 

Theta : Flow angle at the characteristic point (degrees). 

2 2 o 

Entropy : Entropy level at the characteristic point (ft /sec - R). 

Shock Angle : Shock angle at the characteristic point (only prints when 
shock point) (degrees). 

Mach Angle : Mach angle at the characteristic point (degrees). 
Pressure : Static pressure at the characteristic point (psfa). 

Density ; Static density at the characteristic point (slugs/ft ). 
Temperature : Static temperature at the characteristic point (°P). 
Velocity : Velocity of flow at the characteristic point (ft/sec). 

O/F Ratio : O/f ratio at the characteristic point. 

Local Mw : Molecular weight at the characteristic point. 

Local Gamma : Ratio of specific heats at the characteristic point. 

To': Normal shock stagnation temperature at the characteristic point 

(°R). 
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Po : Pitot pressure - the normal shock stagnation pressure at the 

characteristic point (psfa). 

* 

S : Normal shock stagnation entropy at the characteristic point 
(ft^/sec^-°R). 

This is a comparison to the mass flow through the throat. The percent 
change should be near zero; any change is an indication of accumulated 
error in the calculation. For a plume this is a comparison to the mass 
flow at the nozzle exit. The percent change will approach 100% as the 
characteristic lines start further away from the exit plane. The mass 
flow printed out after the first characteristic line is in lb/sec if the 
units of the flow field are in feet. If the flow field is in inches then 
divide the printed out value by 144 to get lb/sec. 

This is a calculation of the components of net momentum thrust and 
pressure thrust at the particular wall point in the nozzle. 

FORCEX: Net force in the axial direction (Ib^) 

FORCEY: Net force in the radial direction. For an axisymmetric 
nozzle FORCEY will be zero (lb^) 

TORQ.Z: The net moment produced by the axial and radial forces 
about the origin of the plume (ft-lb^) 

DELFX: Incremental force in axial direction (lb~) 

DELFY: Incremental force in radial direction (lb^ 

ISPVAC: The vacuum specific impulse (lb^-sec/lb^) 

CF-VAC: The vacuum thrust coefficient. 

ERRORS: and intermediate statement (see Subsection 3.4.2). 
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3.5.2 Statements from Subroutine Errors (Output) 

1. Illegal crossing of right- running characteristics; run is continued 
with point deleted. 

This statement occurs when the newly calculated charac- 
teristic point falls below the right-running characteristic 
line formed by the two base points. 



These points are averaged and the 
average is used to replace the crossed 
points. 


2. Illegal crossing of left-running characteristics; run is continued 

with remainder of line deleted. 

This statement occurs when the new left-running charac- 
teristic line crosses the previous left-running characteristic 
line. 



3-128 


LOCKHEED - HUNTSVILLE RESEARCH & ENGINEERING CENTER 



LMSC/HREC D 162220 -IV 


o 


3. Previously noted errors have propaged to lower boundary or problem 
limits have been reached. Case terminated. 

The program has termirated properly, a shock has inter- 
sected the axis of symmetry, the problem limits set by 
the user have been reached or an error in the mesh con- 
struction which occurred somewhere in the field has 
propagated to the lower boundary. 


4. Lower boundary solution will not converge. 

Program is unable to obtain a solution at lower boundary 
within a preset number of iterations. 


5. Interior solution will not converge. 

Program is unable to obtain a solution at an interior 
point within a preset number of iterations. 


6. Upper boundary solution will not converge. 

Program is unable to obtain a solution at the upper 
boundary within a preset number of iterations. 


7. Shock solution will not converge. Line terminated. 

The program is unable to satisfy both oblique shock 
relations and flow field properties simultaneously within 
a specified number of iterations. 


8. ITSUB WNC in RGMOFP 

Real gas solution of Mach number as a function of 
pressure will not converge in the preset number of 
iterations. 


9. ITSUB WNC in RGVOFM 

Real gas solution of velocity as a function of Mach 
number will not converge in the preset number of 
iterations. 
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10. ITSUB WNC in THETPM 

Unable to balance the last Prandtl-Meyer point 
pressure with back pressure at free boundary or 
flow angle at last point with flow angle at solid 
boundary within the preset number of iterations. 


11. ITSUB WNC in AOASTR 

JU 

Unable to balance the mass flow at A/ A* with 
mass flow at throat within the preset number of 
iterations. 

12. MOCSOL WNC 

The characteristic solution for a new point is 
unable to be reached within the preset number 
of iterations. 


13. SHOCK has problem 


This statement indicates that subroutine (SHOCK) 
is unable to reach a solution. (SHOCK) will reinitiate 
the solution and try again in most cases. When this 
statement is printed out there will also be a 2, 3, 4, 5, 

6 or 7 written out from (SHOCK) which will indicate 
the type of problem. 


2 the shock solution will not converge 

3 the shock has intersected the axis 


4 and 5 not presently used 

6 insufficient downstream information - 
call again with appropriate option 

7 insufficient downstream information - 
call again with appropriate option 


14. ITSUB WNC in TURN 

Iteration for shock angle based on known turning angle 
cannot obtain convergence within specified number of 
iterations. 


15. ITSUB WNC in OVEREX 

Iteration for obtaining shock angle by balancing 
downstream pressure to back pressure. is unable 
to converge within preset number of iterations. 
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16. The following case cannot be found on Master Tape 

The program is unable to find the desired gas case 
among the cases present on the Master Tape. 

17. ITSUB WNC in HYPER 

Unable to balance the pressure at the expansion 
corner with the pressure at the pressure boundary 
within the present number of iterations. 

18. ESHOCK WNC 

Iterative solution of conservation equations across 
a shock would not converge within the preset 
number of iterations. 

19. ITSUB WNC in MASCON 

Iteration for Mach number distribution using con- 
servation of mass flow would not converge within 
the preset number of iterations. 


20. Subsonic Mach number encountered. 

The method of characteristics solution is unable 
to handle subsonic Mach numbers because of the 

use of the te rm^M^ - 1. When this error is 
encountered it is usually the result of the propa- 
gation of some upstream error. 

21. Negative velocity encountered. 

The initial guess for an iterative solution was 
bad which resulted in a negative velocity. May 
also have resulted from some error in the 
boundary conditions. 


22. Characteristic lines diverge, line terminated. 

Unable to find a new characteristic point because 
of a very small difference in angle between two 
characteristic lines. (See figure on following page.) 
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Axis of Symmetry 
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.6 SAMPLE PROBLEM 


Below is a listing of the sample problem. 


Card 1 

5-lV TFSTS RL-IO FNfUNE PB/PC=60.F-5 



Card 2 

1 in 

21 into 

1 o 

1 0113 

o 

o 

1 1 

Card 3a 

0.4*5 r 

>.7 C 

1.0 0 

'.0 




Card 3b 

0.1 

0.325 

0.2 

0.375 

0.3 


0.445 

Card 3c 

0.4 

0.515 

0.5 

0.58 

0.75 


0.725 

Card 3d 

l.o 

0.86 

1 .25 

0.980 

1 .50 


1 .085 

Card 3* 

1 .75 

1.182 

2 *0 

1.28 

2.25 


1.365 

Card 31 

2.5 

1.445 

3.0 

l *59 

3.50 


1.712 

Card 3g 

4.0 

1.83 

4.29 

1.89 




Card 3h 

3 0 30.24 






Card 3i 

2 0 







Card 4 

5- IV TEST* 

-D2/H2-PREL MKS l 

2 



Card 5 

*.o 







Card 6a 

o.o 

1 ? 






Card 7a 

• OOO 

1 1 .558 

1.14376 

3330.81 

23.815995 


Card 7b 

1 .onn 

l 1 .691 

1 • 1 4578 

3141 .66 

13.698737 


Card 7c 

1 .759 

l 1 .299 

U 15768 

2784.51 

4.763100 


Card 7d 

2.148 

l 1 .995 

1 • 1 7286 

2544.81 

2.381*99 


Card 7e 

3.330 

12.094 

1.23567 

1732.2** 

•238160 


Card 7f 

4.237 

12.096 

1.26618 

1254.56 

.04763? 


Card 7g 

4.656 

12.096 

1.28029 

1084.45 

.023816 


Card 7h 

•5.740 

12.096 

1.31156 

760.09 

.004763 


Card 7i 

6.363 

1 2 • 096 

1.72308 

634. 19 

.002382 


Card 7j 

7.888 

12.096 

1.34329 

426.53 

.000476 


Card 7k 

8.71o 

12.096 

1.34937 

353.83 

.000238 


Card 71 

10.834 

12*096 

1.35837 

232.94 

•000048 


Card 6b 

1 .46 

1 2 






Card 7m 

.00" 

10.668 

1.0964 1 

2466*16+6805 

-02 


Card 7n 

1 .Onr% 

10,831 

1.09525 

2389.41+7008 

-02 


Card 7o 

1 .776 

1 1 .14*3 

1 .09457 

2241.21+1361 

-02 


Card 7p 

2.157 

I 1 .345 

1 .09587 

2147.64+6805 

-03 


Card 7q 

3.190 

1 1 .906 

1.12651 

1R08. 12+6805 

-04 


Card 7r 

3.853 

12.086 

1 .22804 

1A52. 72+1 761 

-04 


Card 7s 

4.214 

12.095 

1.26146 

1 267.39+6«05 

-05 


Card 7t 

5.247 

12.096 

1 .29878 

890.26+1361 

-05 


Card 7u 

«.754 

1 2 . 096 

1.31 166 

759.05+65 

-06 


Card 7v 

7.119 

12.096 

1.7352? 

516.42+1361 

-06 


Card 7w 

7.866 

12.096 

1.34310 

428.77+65 

-07 


Card Yx 

0.770 

\ 2.096 

1 . 35478 

283.74+176! 

-07 


Card 8 

0.0 

a 42A 

1 .Ol o 

• n 

0.284 


5.0 

Card 10 

too. n 

• ^ o 

• o r 

.0 20. 

90. 

Card 11 

2.o 0 

• 0 3 

• o r 

• o 

10.0 

0 

'•0 


O 0 


1000 * 

tooo. 


.2 


20 


0 ! 
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Card 1 

Cols . 

Value 

Descriplion 

1-72 


Header card with title of problem. 

Card 2 

5 

1 

ICON(l) - Gas properties will be 
input f-om cards. 

9 

1 

ICON (2) - Right running start line. 

10 

0 

ICON(2) - Straight start line with 
the Mach number given. 

14, 15 

21 

ICON(3) - 21 Start line points 

17-20 

1019 

ICON (4) - The upper boundary will 
be curve fit by the program using 17 
points which describe the nozzle, 
there will be 19 upper wall systems: 
One for the nozzle throat, 17 for the 
nozzle wall and a free boundary 
equation. 

25 

1 

ICON (5) - There is one lower wall 
equation. 

27-30 

0 

ICON(6) - Any shock intersecting 
the lower wall will not be reflected 
automatically. 

35 

1 

ICON(7) - The problem is an 
axisymmetric case. 

37, 38 

01 

ICON(8) - After the solution reaches 
the nozzle lip the output wiT be full 
one line data. 

39,40 

13 

ICON(8) - The output f r the nozzle 
solution up to the lip will be limited, 
three line data. 

41-45 

0 

ICON (9) - The start line does not 
contain a shock point. 

50 

0 

Not presently used. 

54-55 

11 

ICON(ll) - The case number for the 
problem is 11. 

60 

0 

ICON(12) - The first shock encountered 
during the solution will be calculated 

65 

0 

ICON(13) - The viscous boundary layer 


option for the start line v.ill not be used. 
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Card 2 (Cont'd) 


Cols. 

Value 

Description 

69.70 

20 

ICON(l4) - Any Prandtl -Meyer 
expansion will be divided into 
20 "rays." 

75 

1 

ICON(15) - Mesh control will be 
utilized. 

80 

0 

ICON(l6) - No "debugging" printout 
during the shock and pitot pressure 
calculations . 

Card 3a 




1-10 

11-21 

31-40 


31-40 


NOTE: 3a - 3g are used because the nozzle wall curve fit 
option has been selected for the upper boundary- of 
the problem. 

.45 RC - The radius of curvature of the 

nozzle throat is .45 units. 

.3 RT - The throat radius is .3 units. 

0.0 THETA - The throat divergence 

angle is 0.0. The throat equation 
will only be considered up to the 
minimum area of the nozzle. 

0.0 XO - The nozzle throat equation will 

not be shifted. 


Card 3b 


1-10 

.1 

11-20 

.325 

21-30 

.2 

31-40 

.375 

41-50 

.3 

51-60 

.445 


XIN(l) - Axial location of first point 
describing nozzle wall. 

YIN(l) - Radial location of first point 
describing nozzle wall. 

XIN(2) 

YIN(2) 

XIN(3) 

YIN(3) 


Cards 3c-3g 


These cards contain the coordinates of the rest of the 17 points de- 
scribing the nozzle wall. 
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Card 3h 



Cols. 

Value 

Description 

1 

3 

rWALL - The last upper boundary 
equation is a free boundary equation. 

5 

0 

ITRANS - No discontinuity follows 
this boundary. 

11-20 

30.24 

PINF - The nozzle will expand to a 
free boundary back pressure of 30.24 
psfa. 

61-70 

1000. 

XMAX - The maximum axial coordi- 
nate for which this equation applies 
is 1000. units. 

Card 3i 



1 

2 

IWALL - The lower wall equation is 
a polynomial equation. 

5 

0 

ITRANS - No discontinuity follows 
this equation. 

11-60 

0.0 

WALLCO(l)-(5) - The lower boundary 
equation is a straight line with zero 
slope passing through the origin of 
the system. 

61-70 

1000. 

XMAX - The maximum axial coordi- 
nate for which this equation applies 

Card 4 


is 100 units. 

1-24 

S-IV Tests 
-02/H2 PR EL 

ALPHA - Gas name and description. 

30-32 

MKS 

UNITS - The gas data is input in 
metric units. 

40 

1 

IOF - There is 1 O/F table of gas data. 

45 

2 

IS - There are two entropy tables con- 
taining gas data. 

Card 5 



1-10 

5.0 

OFRAT - The O/F ratio for this data 
is 5.0. 

Card 6a 



1-10 

0.0 

STAB(l) - The entropy level of the 
first entropy table. 
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Card 6a (Cont'd ) 


Cols. 

Value 

Description 

19-20 

12 

IVTAB - There are 12 velocity or 
Mach number "cuts" of data in this 
entropy table. 

Card(s) 7a 



1-10 

0.0 

TAB (I, J,K,1) - Mach number for the 
first Mach number cut in the entropy 
table. 

11-20 

11.558 

TAB(I,J,K,2) - Molecular weight of 
the gas at this Mach number. If the 
units on Card 4 were ENG then this 
should be the gas constant R. 

21-30 

1.14376 

TAB (I, J,K,3) - GAMMA for the gas 
at this Mach number. 

31-40 

3330.81 

TAB(I, J,K,4) - Static temperature at 
this Mach number, °K. 

41-50 

23.815995 

TAB(I,J,K,5) - Static pressure at this 
Mach number, Atmospheres. 

Cards 7b- IS. 



These cards contain similar information to that cm Card 7a except 

corresponding to different Mach number "cuts." 

Card 6b 



1-10 

1.46 

STAB (2) - The entropy level of the 
second entropy table (cal/gram-°K) 

19-20 

12 

IVTAB - There are 12 velocity or 
Mach number "cuts" of data in this 
entropy table. 

Cards 7m-7x 



These cards 

contain similar data 

as Cards 7a-7 i except that this data 

was generated for 

an entropy level corresponding to 1.4 cal/gram-°K. 

Card 8 



1-10 

0.0 

CORLIP(l) - The axial coordinate of 


the upper limit of the start line is 0.0. 
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Card 8 (Con'd) 


Cols . 

Value 

Description 

11-20 

.0428 

CORLIP(2) - The axial coordinate of 
the lower limit of the start line is 
.0428 units. 

21-30 

1.01 

CORLIP(3) - The Mach number of the 
flow along the start line. 

31-40 

0.0 

CORLIP(4) - There is zero entropy 
level at the start line. 

41-50 

.284 

STEP(5) - The area of the nozzle 
throat. 

51-60 

5.0 

CORLIP(5) - The O/F ratio along the 
start line. 

61-70 

0.0 

STEP(6) - The program will set the 
minimum AP for discontinuing shock 
calculations at .001. 

Card 10 

1-10 

100. 

CUTDAT(l) - The maximum radial 
coordinate considered for this prob- 
lem will be 100. units at the upstream 
cutoff point. 

11-20 

0.0 

CUTDAT(2) - Each characteristic line 
will cut off for any axial value less 
than 0.0. 

21-30 

0.0 

CUTDAT(3) - The angle the upper 
cutoff makes with the horizontal is 0.0 

31-40 

0.0 

CUTDAT(4) - The minimum radial 
coordinate considered for this prob- 
lem will be 0.0 at the downstream 
cutoff point. 

41-50 

20.0 

CUTDAT(5) - The maximum axial 
coordinate for which calculations 
will be made. 

51-60 

90.0 

CUTDAT(6) - The angle the down- 
stream cutoff line makes with the 
horizontal is 90.0 degrees. 

61-70 

.2 

STEP(3) - Anytime the axial distance 


between two characteristic lines ex- 
ceeds .2 units a new left ruraing char- 
acteristic line will be st-rted’. 
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Card 11 


Col(s) 

Value 

Description 

1-10 

2.0 

DIV(l) - Maximum change in absolute 
length between any two adjacent char- 
acteristic points is 2.0 units. 

11-20 

0.0 

DIV(2) - No points will be deleted be- 
cause of the minimum change in abso- 
lute length. 

21-30 

3.0 

DIV(3) - Maximum change in Mach 
number between any two consecutive 
characteristic points is 3.0. 

31-40 

0.0 

DIV(4) - No points will be deleted 
because of the minimum change in 
Mach number. 

41-50 

10.0 

DIV(5) - Maximum change in flow 
angle between any two consecutive 
characteristic points is 10°. 

51-60 

0.0 

DIV(6) - No points will be deleted 
because of the minimum change in 
flow angle. 
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Appendix 

SOLUTION TECHNIQUES 



LOCKHEED - HUNTSVILLE RESEARCH & ENGINEERING CENTER 



LMSC/HREC D1 62220-1 V 


Appendix 

MESH CONSTRUCTION FOR INTERNAL FLOW 

The calculations described previously are point or small region 
solutions. Some process must be defined which successively employs 
the proper calculation at the proper time in order to describe the entire 
field. In order to facilitate a description of the mesh construction process 
let $ represent the total knowledge of flow properties at a point in the field. 
Also let the expression 


* = * 1 * 1 -*2 * m ) 


stand for properties at a new point which are computed as a function ifj/j of 
(m) other points. There will be basically six such functional operations 
xlt ^ jpg, <J/ S tf^M' s ^ an( ^ f° r input point, interior point, bound- 


ary point, attached shock point, shock, and Prandt-Meyer points. In addi- 
tion the superscript (u) will indicate that the operation is to be performed 
in the presence of an upper boundary while (L) indicates a lower boundary. 


Due to the complexity of handling multiple shock waves, a single 
shock wave restriction will be imposed. This shock wave is arbitrarily 
chosen to be of the right running family. This type of problem will most 
frequently occur in cases of interest. If a left running shock wave occurs 
the problem is simply inverted. 

The choice of right running shock waves also dictates that left running 
characteristic lines be followed in the calculation. This allows one to re- 
tain a minimum of information, i. e., a known characteristic line (hereafter 
referred to as (j)) and a line in the process of being computed (k). 
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To begin the problem all necessary boundary conditions must, of 
course, be supplied. In addition a starting line containing N points which 
are designated ^ (n = 1, . . . N) must be supplied. 

Figure A-l illustrates a flow field in which there are no discontinuities 
and in which the mesh construction is terminated when the region of interest 
has been computed. 



In region I the left running characteristic lines initiate as input points 
and the mesh construction may be described by; 


■ 

w 

i = 1 


♦i.k = < 


i = 2 . . . (2n-2) 

(A 


( ^B^i-l.k’ *i-2,j> 

i = 2n-l 
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where n varies from 1 to N and 4>. , represents the flow properties at the 
th * 

i point on the k line . For instance , in calculating the fourth point on the 
fourth line shown in the figure, line three is known in its entirety and line 
four up to and including the third point is known. The above set of relations 
says that point three on line four and point three on line three will deter- 
mine, through the interior point solution, the next point (four) on line four. 


For region II we have; 


i = 1 


, 

x, k 


*5 «*!,!■ *2,j> 

1 = 2 -”-< 2N - 2 > 


(A-2) 


1 = 2N -‘ 


As a new line becomes completely defined it may be referred to as j 
and the process continued indefinitely. 


It is possible to combine regions (1) and (2) into a more general 
scheme if a variable ij^ is defined which takes on the value (1) in region I 
and (0) in region II. At this time the number of points on the j line (i^, ) 

j 

and the number of points on the new line (i T ) are defined. Then; 


4 >. 


i, k 


N 


W + » 2 .j> 

V 


^ik.k’ *ij, j ) iji = i-2i +1 


i = 1 


i — 2, . . . i_ - 1 
X k 


(A-3> 


if/ 1 (* & \ ik = i-1 

o * xk, k’ ij, j 1 ij = i-2i, 


N 


i = x. 


but 


x <= i + Zi 
k J 
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Obviously i^ would have been initialized to (-1) prior to the start oi the 

j 

calculation. When the line is fimsned i_ is set to the current value of i^ . 

J k 

Thus the process for computing the entire flow field for such a simpli- 
fied case is described by the set of expressions (A- 3). In general, however, 
discontinuities will arise so that a more flexible description is necessary. 

If, by some process, points were discarded from the (j) and (k) arrays and 
the number of points lost is ig and ig respectively then (A- 3) becomes; 

j k 


i-.tf' {$ ) + (l-i_ T )i/£ (<f> <l> .) 

N o' n N B 1 , j 2, j 


i = 1 


ik = i- 1 


^i,k = ^ ^i^ik, k* ^ij, ij = i-2i N +l+i^+i 




i = 2 i - 1 (A- 


^B^ik.k’ *ij,j* ij = i-2i N + ig +i g V 1 " ln 

j k 


but 


i T. + 2i N- i 8- 
J J k 


where the tilda over ig^ and ig^ indicates that current values are to be used. 
These variables are reset to zero at the beginning of each new line. 


To illustrate this, imagine that points (1) and (2) have been computed 
on line 11 and that after point(3 / ) had been computed in the normal fashion 
it was necessary to discard it. The next point to be computed would then 
be (3*) but if for some reason it was necessary to discard point f» or. the 
(j) line then the point ( 3 ^) would not exist. Therefore a point has been 

deleted on each line (i^ = i^ =1) and the diagram and the set of equations 

j k 

(4) indicate that point (6) on line (j) and point (2) on line (k) would be used 
in the computation of (3) on line (k). 
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It is now possible to include a shock wave into the logic scheme. 
Since it i6 a mathematical requirement that characteristic lines of the same 
family as the shock are continuously intercepted by it, the ability to discard 
points was necessary. If this was the only mechanism for discarding points 
then the logic process would be; 


4*. , 

x, k = 


« =i-2i N + 1 | 




1 ~ *»k " 1J k 
1 =i »k- i Sk + 1 


i = I 


i = 2, . . . i 


(A- 5) 


^ik. k’ ^ij, ^ V’l. ..... li-i -i* 

1J = 1_2i n u \ ( k 5 


+ 2, . . .i_ -1 
k X k 


if} 1 (A A | 

B 1 ik, k* ij, j ; ij = i-2i^+ ig + ig 

k j 


ik = i-1 


i = i. 


where 


i T, ^ i T. + ‘ i 8- “ % 
k J J k 


and where i is defined in much the same fashion as i^, , which is; 

S. * i 

k k 


+ 2i 




where i 

s . 
J 


is the location of the upstream shock point on the (j) line. 
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Figure A-2 illustrates the mesh construction when a shock wave 
is present. 



In this example i. T = 0 and i = 8 so that i =7. Also i = 16 so 

JN S . S» X . 

j k j 

that i_ would normally also be 16. The (k) line is computed up to point 7 
A k 


and the shock solution is then employed. In this case the shock solution 
finds that three points of the (k) line fall downstream of the shock (minimum 
is one) while two right running lines (points 10 and 11 on the (j) line) also 
are intercepted by the shock wave. Thus in this example ig =2 and ig = 2. 

k j 


The set of equations (A-5) then says that the double shock point should be 

points 5 and 6 on the (k) line and that the total number of points on the (k) 

line has decreased to 12. Note also that the value of i to begin the next 

s * 

line must be change to 5. ^ 


So far no mention of how the shock wave begins has been made. 
There are two types of shock waves considered; the attached shock wave 
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which arises due to the flow being forced to negotiate a compression corner 
on the upper boundary, and the envelope shock. The first of these is easily 
detected from the boundary conditions and is initially of finite strength. The 
second type is detected by a mathematical discontinuity in the mesh con- 
struction (crossing of right running lines) and is initially of zero strength 
i.e., a Mach wave. An example of the compression corner solution is given 
in Figure A- 3. 



The computation of the (k) line is completed without any prior know- 
ledge that a compression corner exists. A check is made after the bound- 
ary solution and the boundary information indicates that a compression 
corner must be treated. A linear interpolation is performed between the 
boundary point on line (j) and the fictitious boundary point on line (k) in 
order to determine the flow properties at point u. An oblique shock calcu- 
lation is made where the turning angle is known. Using this point and point 

6 (i_ -1) a new virtual point (/) is computed, i is set to i_ and the shock 
k S k k 

solution illustrated in Fig. A-4 is employed. This shock solution completes 

the (k) line in the proper fashion. In this example i^ = 6 and the next line 

is computed as previously discussed. ^ 
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Figure A-4 - Solid or Free Upper Wall Interacting with Reflected Shock 
Wave or Attached Shock Wave with Insufficient Downstream 
Information 

The envelope shock is detected by a crossing of right running char- 
acteiistic lines as shown in the figure below. 



Figure A-5 

In this example point (5) on the (k) line is found to fall in a previously 
described region (the region between points (3) and (4) on the (k) line). This 
discontinuity in the solution is interpreted as a shock wave. If the grid size 
were chosen small enough the shock wave would initially be of zero strength. 
Point (5) on line (j) is chosen to be a point which lies on the shock wave and the 
shock solution is employed. The results of this solution are stored in the normal 
fashion on the (k) line. Obviously the only difference between this situation and 
treatment of a previously developed shock wave is to modify the (j) line such 
that it appears to the logic scheme as though a shock wave crossed the (j) 1 ne 
at point (5). 
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Figure A-6 illustrates the mesh construction in the vicinity of an expansion 
corner. 



Figure A-6 


In this case point (5) on line (k) is expected to be a boundary point. 

It is discovered however, that an expansion corner must be negotiated by 
the (k) line. A point (6) on the (j) line is found by interpolation. A Prandtl- 
Meyer calculation is employed and the fan of points is stored in the (j) line 
above (6). The total number of points to be expected on the (k) line is in- 
creased accordingly and the normal logic scheme will now complete the 
new line. The next line is calculated in the standard fashion. 

An expansion corner on a lower wall is somewhat a more complicated 
situation. Since the calculation no longer utilizes the input line, the lower 
wall expansion fan may be stored in this area. The set of relations, Equa- 
tion (A-5) is modified to that of Equation (A-6). 
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I *2,j> 


ik = i-1 


^L^ik.k’ *ij, j* ij = i + (i^-2; i^. 4- 1 | 


% ( *k’ *j> 


<t>.. 

!J = 


1 =1 -lc 

s k S k 

i = i -irs + 1 

s k °k 


i = 2, 


|^*ik,k* *ij, j* ij = i + (i J -2)i N + ig + i g 1 i=i 

k j ) 


*1 (*;u u.*:: ;> « I \~L 


B ik, k’ ij, j ij = i + (i | -2)i N + i 


V‘«- I i = i ' 

k J 


where 


*T, 1 T. + (2_1 i ) i N _1 5-" 1 8, 

k j J k 


and where 


‘s. <= i s. + < 2 -V i N - 1 

k J 


for which i^ = i = 1 until all fan points are used up. 
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Mesh Control 

The mesh control utilized by the program is a means of controlling 
the minimum and maximum changes in distance, Mach number and flow 
angle between any two consecutive characteristic points on both left and 
right running characteristic lines. The advantage of mesh control occurs 
in plumes and nozzles where a great deal of expansion is encountered. The 
characteristic mesh size expands to a degree where large spaces may be 
encountered in the plume and insufficient data is present. 

Mesh control is employed by checking on the absolute difference in 
Mach number and flow angle and the distance between two points. These 
checks are made between both the left and riglit running base points and the 
newly calculated characteristic point. If these differences do not satisfy 
the input criteria the new point will either be deleted or another pointed 
added. The fixed point equations and counting schemes within the program 
are altered to account for points which are added and deleted. 

The following two figures are examples of insertion and deletion of 
points as a result of mesh control criteria. 


A-l 1 


LOCKHEED - HUNTSVILLE RESEARCH & ENGINEERING CENTER 



LMSC/HREC D162220 -IV 


Example of Insertion 



Figure A-7 


If (ABSfQ^-S^)) is greater than DIV(5) 
or (ABSfM^-M^)) is greater than DIV(3) 
or (ABS(AL <2 ^)) is greater than DIV(l) 
a new left running point 3 A is added 
which starts a new left running charac- 
teristic line. The following is an 
example of the message printed out for 
this case: 

EXTRA LRC NO. 1 HAS BEEN ADDED 
BETWEEN LINES 43 AND 44 BECAUSE 
OF A CHANGE OF . 20468+00 DID NOT 
SATISFY THE DISTANCE MESH CON- 
STRUCTION BETWEEN POINTS 4 AND 3. 


If (ABSfO^-O^)) is greater than DIV(5) 
or (ABStM^-M^)) is greater than DIV(3) 
or (ABS(AL.^ ^)) is greater than DIV(l) 

a new right running point 3 will be added between point 2 and old point 3. Old 
point 3 then becomes point 4. The following is an example of the message printed 
out for this case: 

A NEW RRC HAS BEEN ADDED BETWEEN POINTS 2 AND 3 ON LINE 
44 BECAUSE A CHANGE OF .20140+00 DID NOT SATISFY THE FLOW 
ANGLE MESH CONSTRUCTION. 


It is possible for more than one left running characteristic to be added as shown 
when a new left running point 6B was added between points 5A and 6 on line 44A 
and 44. 
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Example of Deletion 



If (ABSfO^-Q^) * s * ess than DIV(6) 
or ABS(0 2 -0 1 ) is less than DIV(6) 
or ABSfM^-M^ is less than DIV(4) 
or ABS(M 2 -Mj) is less than D1V(4) 
or ABS(ALij 2 ) i® less than DIV(2) 
or ABS(AT , is less than DIV(2) 

point (2) is deleted and a new point 2 is calculated using point 4 
on line 40 and point 1 on line 41. 

The following message will occur for this example: 

THE POINT 2 ON LINE 41 HAS BEEN DELETED BECAUSE A 
CHANGE OF 1.5 DID NOT SATISFY THE MACH NUMBER MESH 
CONSTRUCTION. 
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Automatic Shock Reflection 

The method of automatically reflecting a shock begins when the shock 
intersects the axis of symmetry. The downstream shock point and the last 
point on a specified boundary equation are connected with a straight line. 
The flowfield is then searched for flowfield properties which fall at certain 
points on the line as determined from the number of specified number 
starting line points. If no point is located which falls on the line, an inter- 
polation is performed to determine the flowfield properties at the point. 
Once the specified number of points are found the properties at each point 
are stored in the starting line array. The problem and boundary equation 
are then inverted and the problem in reinitiated. 
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Viscous Start Line Option 

The viscous start line option is for use in setting up a supersonic viscous 
one -dimensional start line at a nozzle exit plane. Subroutine VISCOUS is called 
with N, the inverse of the power of the profile, NBL, the number of points to 
be calculated in the viscous portion of the start line, XL., the characteristic 
length (ft. - usually the nozzle wall length) and CU, a conversion factor for 
converting the characteristic length XL into the units of the flow field. 

A flat plate boundary layer formulation was used to determine the boundary 
layer properties. Since the local flow properties vary greatly through the 
boundary layer of supersonic flow an average value of pertinent flow properties 
should be determined. A reference temperature obtained by Equation 1 (Ref. 1) 
is used to determine 

T ref = T ex * - 5 <' 9T ° ' T e*> + - 22 < T o ' T «> (1) 

an average viscosity, where the viscosity is obtained from Equation 2. 

= (2.27 (T ref X * 5 ) 1. ” 8 )/(T ref + 198.6) (2) 

Using these reference conditions a local Reynolds number per foot is obtained 
using Equation 3. 

R = P V XLAa RT . ,-ix 

e ex ex ref (3) 

The Reynolds number determined from Equation 3 is checked to determine 
if the boundary layer is laminar or turbulent. A Reynolds number less than 
1 million is laminar and a Reynolds number greater them 1 million is turbulent. 

Once the flow regime has-been determined it is then necessary to determine 
the boundary layer thickness. The boundary layer thickness is taken as that 
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distance from the surface at which the local velocity reaches 99 percent of the 
free stream velocity. Equations 4 and 5 are used to determine the boundary 
layer thicknesses of a laminar or turbulent boundary layer, respectively. 


S L = 5.2 XL*CU/R e 5 (1 + .732FL) 

(4) 

& T , = .0188 XL*CU*M*‘ 25 (1 + .176F)/R; 14 

(5) 


FLand F are terms used to adjust the boundary layer thickness depending on 
whether or not the nozzle is two-dimensional or axisymmetric. 

After the boundary layer thickness has been determined the subsonic 
portion of the boundary layer is eliminated from the calculated thickness. 

The starting line points (NBL) are then distributed through the boundary layer 
using an empirical equation of the form of Equation 6 (Ref. 2). 



N is assumed to be a 2 if the boundary layer is laminar. If the boundary layer 
was found to be turbulent the value of N input by the programmer is used, 

(N = 7 to 9). 

After the boundary layer points are distributed through the boundary layer 
they are transferred to a right running characteristic and the rest of the start 
line points are adjusted to the edge of the boundary layer. 
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Symbols for Viscous Start Line Option 


CU 

FL 

F 

M 

e 

N 
P 
R 
R 


ex 


ref 

v 

V ex 

XL 

y 

5 X 
* 


el 


A*- 


Conversion factor 

1 for axisymmetric nozzle, 0 for two-dimensional nozzle 

0 for axisymmetric nozzle, 1 for two-dimensional nozzle 

Mach number at edge of boundary layer 

Reciprocal of the exponent of the profile 

Static pressure through the boundary layer 

Gas constant 

Local Reynolds number 

Total temperature at edge of boundary layer 
Eckert's reference temperature 
local velocity 

Velocity at edge of boundary layer 
Nozzle wall length - ft. 

Distance from nozzle wall 
Laminar boundary layer thickness 
Turbulent boundary layer thickness 

Boundary layer thickness with subsonic portion eliminated 
Reference viscosity 


References for Viscous Start Line Option 


1. Eckert, E. Ri, "Survey of Heat Transfer at High Speeds," WADC 54-70, 
Wright Air Development Center, April 1954. 

2. Schlechting, H., "Boundary Layer Theory," McGraw-Hill Book Company, 
Inc., New York, New York, I960. 


A-17 


LOCKHEED -HUNTSVILLE RESEARCH & ENGINEERS. .ENTER 



